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Nomenclature 


AR 

=  Aspect  ratio 

rrc 

=  Inboard  dihedral  angle,  deg 

S 

=  Wing  wetted  area,  m2 

Pet 

=  Outboard  dihedral  angle,  deg 

CR 

=  Root  chord  length 

=  Freestream  Mach  number 

b 

=  Semi  span  length 

Re 

=  Reynolds  number 

Arc 

=  Inboard  taper  ratio 

a 

=  Angle  of  attack,  deg 

Act 

=  Outboard  taper  ratio 

W 

=  Yaw  angle,  deg 

Arc 

=  Inboard  sweep  angle,  deg 

CL 

=  Lift  coefficient 

Act 

=  Outboard  sweep  angle,  deg 

cD 

=  Drag  coefficient 

BP  Inboard 

=  Inboard  break  point 

Cdo 

=  Drag  coefficient  at  zero  lift 

BP  Outboard 

=  Outboard  break  point 

L/D 

=  Lift  to  drag  ratio 

Cf 

=  Friction  coefficient 

SrTr 

=  Spar  Root  Thickness  Taper  Ratio 

cranki 

=  Crank  Location 

WsRt 

=  Wing  Skin  Root  Thickness,  m 

Ns 

=  Number  of  Spars 

WstTr 

=  Wing  Skin  Thickness  Tip  Taper  Ratio 

Nr 

=  Number  of  Ribs 

WsTre 

=  Wing  Skin  Thickness  Edge  Taper  Ratio 

Rrt 

=  Rib  Root  Thickness,  m 

Sc 

=  Spar  Cap  Root  Area,  m2 

RrTr 

Srt 

=  Rib  Root  Thickness  Taper  Ratio 
=  Spar  Root  Thickness,  m 

Re 

=  Rib  Cap  Root  Area,  m2 
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Summary 


The  overall  objective  of  this  project  was  to  realise  the  practical  application  of  Hierarchical  Asynchronous 
Parallel  Evolutionary  Algorithms  for  Multi-objective  and  Multidisciplinary  Design  Optimisation  (MDO)  of 
UAV  Systems  using  high  fidelity  analysis  tools.  The  project  looked  at  the  aerodynamics  and  structure  of 
two  production  UAV  wings  and  attempted  to  optimise  these  wings  in  isolation  to  the  rest  of  the  vehicle. 
The  project  was  sponsored  by  the  Asian  Office  of  the  Air  Force  Office  of  Scientific  Research. 

The  two  vehicles  wings  which  were  optimised  were  the  Northrop  Grumman  Global  Hawk  (GH)  and  the 
General  Atomics  Altair  (Altair).  The  optimisations  for  both  vehicles  were  performed  at  cruise  altitude  with 
MTOW  minus  5%  fuel  and  a  2.5g  load  case.  The  GH  was  assumed  to  use  NASA  LRN  1015  aerofoil  at  the 
root,  crank  and  tip  locations  with  five  spars  and  ten  ribs.  The  Altair  was  assumed  to  use  the  NACA4415 
aerofoil  at  all  three  locations  with  two  internal  spars  and  ten  ribs.  Both  models  used  a  parabolic  variation  of 
spar,  rib  and  wing  skin  thickness  as  a  function  of  span,  and  in  the  case  of  the  wing  skin  thickness,  also 
chord. 

The  work  was  carried  out  by  integrating  the  current  University  of  Sydney  designed  Evolutionary  Optimiser 
(HAPMOEA)  with  Computational  Fluid  Dynamics  (CFD)  and  Finite  Element  Analysis  (FEA)  tools.  The 
variable  values  computed  by  HAPMOEA  were  subjected  to  structural  and  aerodynamic  analysis.  The 
aerodynamic  analysis  computed  the  pressure  loads  using  a  Morino  class  panel  method  code  named 
PANAIR.  These  aerodynamic  results  were  coupled  to  a  FEA  code,  MSC.Nastran8  and  the  strain  and 
displacement  of  the  wings  computed.  The  fitness  of  each  wing  was  computed  from  the  outputs  of  each 
program. 

In  total,  48  design  variables  were  defined  to  describe  both  the  structural  and  aerodynamic  properties  of  the 
wings  subject  to  several  constraints.  These  variables  allowed  for  the  alteration  of  the  three  aerofoil  sections 
describing  the  root,  crank  and  tip  sections.  They  also  described  the  internal  structure  of  the  wings  allowing 
for  variable  flexibility  within  the  wing  box  structure.  These  design  variables  were  manipulated  by  the 
optimiser  such  that  two  fitness  functions  were  minimised.  The  fitness  functions  were  the  overall  mass  of 
the  simulated  wing  box  structure  and  the  inverse  of  the  lift  to  drag  ratio.  Furthermore,  six  penalty  functions 
were  added  to  further  penalise  genetically  inferior  wings  and  force  the  optimiser  to  not  pass  on  their  genetic 
material. 

The  results  indicate  that  given  the  initial  assumptions  made  on  all  the  aerodynamic  and  structural  properties 
of  the  HALE  and  MALE  wings,  a  reduction  in  mass  and  drag  is  possible  through  the  use  of  the 
HAPMOEA  code.  The  code  was  terminated  after  300  evaluations  of  each  hierarchical  level  due  to  plateau 
effects.  These  evolutionary  optimisation  results  could  be  further  refined  through  a  gradient  based  optimiser 
if  required.  Even  though  a  reduced  number  of  evaluations  were  performed,  weight  and  drag  reductions  of 
between  10  and  20  percent  were  easy  to  achieve  and  indicate  that  the  wings  of  both  vehicles  can  be 
optimised. 
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1  Introduction 


The  development  and  use  of  Unmanned  Aerial  Vehicles  (UAVs)  for  military  and  civilian  applications  is 
rapidly  increasing.  Similar  to  their  manned  counterparts,  the  challenge  is  to  develop  trade-off  studies  of 
optimal  configurations  to  produce  high  performance  aircraft  that  satisfy  mission  requirements.  There  are 
difficulties  with  these  concepts  due  to  the  varied  nature  of  the  missions  that  have  to  be  performed.  These 
complexities  in  engineering  design,  and  the  more  demanding  industrial  requirements,  have  pushed  the  need 
to  speed  up  the  development  of  robust  and  fast  numerical  techniques  to  overcome  difficulties  associated 
with  traditional  deterministic  optimisers. 

In  aerospace  engineering  design  and  optimisation  the  designer  is  usually  presented  with  a  problem  which 
involves  not  just  one  objective,  but  numerous  objectives  and  multi-physics.  Hence  a  systematic  approach, 
which  is  regarded  as  Multidisciplinary  Design  Optimisation  (MDO)  and  which  accounts  for  the  coupling 
between  the  variables  and  objectives,  is  required.  Problems  in  aeronautics  usually  involve  a  number  of 
disciplines  and  objectives;  therefore  the  search  space  can  be  multi-modal,  non  convex  or  discontinuous. 
Wing  design  is  an  example  of  a  multi-objective  MDO  problem  where  there  is  a  strong  interaction  between 
the  aerodynamics  and  the  forced  structural  response. 

The  goal  of  the  present  report  is  to  address  this  issue  from  a  multi-objective  and  multidisciplinary  aircraft 
design  optimisation  (MDO)  standpoint. 

In  today's  world  of  advancing  technology,  engineers  are  faced  with  the  problem  of  designing  increasingly 
complicated  multidisciplinary  systems.  This  is  a  difficult  task,  as  these  systems  not  only  involve  coupling 
amongst  the  different  physics  involved,  but  also  a  large  number  of  variables  and  a  series  of  objectives  and 
constraints.  At  the  same  time,  engineers  need  to  optimise  and  address  several  requirements  which  include 
the  reduction  of  the  time  spent  on  design  and  reducing  the  cost  of  research  and  development,  while 
improving  the  performance,  reliability,  quality  and  safety  of  the  product  or  process  under  consideration. 

These  complex  interactions  have  generated  a  growing  interest  in  the  area  of  multi-objective  and 
Multidisciplinary  Design  Optimisation  (MDO).  Here  the  designer  is  interested  not  only  in  a  single  global 
optimal  solution  but  in  a  set  of  solutions  that  represent  the  trade-off  between  the  different  objectives.  MDO 
refers  to  an  approach  to  formalise  a  design  process  which  accounts  for  the  interaction  between  the  different 
physics  involved,  while  optimising  for  a  number  of  objectives  and  constraints. 

When  applied  to  aeronautics,  the  necessity  of  optimisation  and  MDO  is  clear,  given  that  even  a  very  small 
improvement  in  weight  or  a  reduction  in  aerodynamic  drag  will  have  a  tremendous  impact  on  the  overall 
performance  of  the  design.  The  application  of  MDO  to  the  design  and  optimisation  of  aerospace  vehicles 
can  be  represented  as  in  Figure  1.  Cleary  in  MDO  a  number  of  disciplines  (aerodynamics,  structures, 
propulsion,  aero-acoustics,  etc)  are  present  and  interact. 

The  task  of  the  designer  or  team  of  designers  is  to  develop  a  solution  that  conforms  to  all  the  disciplines 
while  satisfying  the  requirements  and  constraints.  When  optimisation  is  intended,  the  different  objectives 
(for  example,  aerodynamic  performance,  purchase  price,  take-off  weight)  need  to  be  considered  in  order  to 
find  an  optimal  solution  or  set  of  solutions. 

A  common  approach  for  optimisation  is  the  use  of  aggregating  functions  in  which  different  weights  are 
assigned  to  each  objective.  The  problem  with  this  approach  is  that  the  weight  for  each  objective  needs  to  be 
known  in  advance.  Another  approach  is  to  compute  or  produce  a  set  of  solutions  in  what  will  be  referred  to 
in  this  thesis  as  a  Pareto  optimal  front  or  surface.  This  Pareto  optimal  front  represents  the  optimal  set  of 
non-dominated  solutions  and  the  trade-off  between  the  objectives  and  disciplines  involved. 
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The  process  of  MDO  involves  the  use  of  several  analysis  tools  such  as  Computational  Fluid  Dynamics 
(CFD)  software  or  Finite  Element  Analysis  (FEA)  and  also  optimisation  tools.  These  analysis  tools  are 
under  constant  development  and  have  reached  a  point  where  a  confident  application  to  aeronautical  design 
in  conjunction  with  MDO  is  possible  (Mason  1998;  Agarwal  1999;  Thomas  1999) 

However,  there  are  limitations  on  the  application  of  these  tools  within  MDO  at  an  industrial  level,  due  to 
the  computational  expense  involved.  A  single  high-fidelity  Navier-Stokes  CFD  computation  around  an 
aircraft  wing,  for  example,  might  take  several  hours  on  a  supercomputer.  Therefore,  the  continuing 
challenge  has  been  to  develop  methodologies  such  as  Design  of  Experiments  (DOE)  (Giunta  1997), 
approximation  methods  and  variable  fidelity  models  that  combine  and  use  different  fidelity  analysis  tools 
during  the  design  and  optimisation  process  to  minimise  the  computational  expense. 

While  the  area  of  traditional  optimisation  tools  for  a  single  discipline  is  quite  mature,  the  area  of  robust 
optimisation  tools  and  approaches  for  MDO  is  still  at  the  initial  stages  of  development  (Alexandrov  and 
Lewis  2000),  (Sobieszczanski-Sobieski  1997;  Bartholomew  1998) 

The  conventional  approach  in  aeronautical  design  and  MDO  has  been  the  use  of  traditional  deterministic 
opt i misers.  These  optimisers  are  efficient  for  finding  optimal  global  solutions  if  the  objective  and 
constraints  are  differentiable.  However,  robust  alternative  numerical  tools  are  required  if  a  broader 
application  of  the  optimiser  is  desired,  or  the  problem  is  multi-modal,  involves  approximations,  is  non- 
differentiable  or  involves  multiple  objectives  and  physics,  as  is  usually  the  case  in  the  design  of  complex 
multidisciplinary  systems  in  aeronautics. 

A  relatively  new  technique  for  optimisation  is  the  use  of  Evolutionary  Algorithms  (EAs).  These  are  based 
on  Darwinian  theories  of  evolution,  where  populations  of  individuals  evolve  over  a  search  space  and  adapt 
to  the  environment  through  the  use  of  different  mechanisms  such  as  mutation,  crossover  and  selection.  EAs 
require  no  derivatives  or  gradients  of  the  objective  function  and  have  the  capability  of  finding  globally 
optimum  solutions  amongst  many  local  optima.  They  are  easily  executed  using  parallel  computing 
techniques  and  can  be  adapted  to  arbitrary  analysis  codes  without  much  major  modification.  Another  major 
advantage  of  EAs  is  that  they  can  tackle  multi-objective  problems  directly.  Together  these  characteristics 
give  EAs  substantial  advantages  over  more  conventional  deterministic  approaches. 
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Interest  in  EAs  for  problems  in  engineering  and  aeronautics  has  grown  substantially  in  the  past  fifteen 
years.  These  methods  have  been  successfully  applied  to  different  aeronautical  design  problems  including 
airframe,  wing,  aerofoil  and  rotor  blade  design,  (Obayashi  1998),  (Oyama,  Liou  et  al.  September,  2002), 
(Parmee  and  Watson  July  1999),  (Raymer  2003;  Gonzalez  2005). 

The  application  of  EAs  for  MDO  problems  has  been  limited.  This  is  mainly  because  one  of  the  drawbacks 
of  EAs  is  that  they  are  slow  when  compared  to  traditional  deterministic  methods,  as  they  require  a  larger 
number  of  function  evaluations  to  converge  to  an  optimal  solution. 

Hence  the  continuing  challenge  in  evolutionary  optimisation  has  been  to  reduce  the  number  of  function 
evaluations  and  the  computational  expense.  To  achieve  this,  several  approaches  have  been  proposed;  these 
include  a  combination  of  variable  fidelity  models,  parallelisation  strategies  and  hybridisation  techniques 
(Coello,  Veldhuizen  et  al.  2002),  (Deb  2003),  (Kim  1 997). 

This  report  describes  the  theory  and  application  of  a  methodology  for  multi-objective  multidisciplinary 
design  of  UAV  Wings.  The  method  is  based  on  a  robust  evolutionary  optimiser,  an  aero-structural  solver, 
parallel  computing,  asynchronous  evaluation  and  a  hierarchical  topology  of  different  fidelity  solvers  that 
reduce  the  computational  expense  for  multi-objective  and  MDO  problems.  The  method  is  applicable  to 
single  and  multi-objective,  inverse  or  direct  complex  engineering  problems  that  can  be  multi-modal, 
involve  approximations,  are  non-differentiable,  with  convex,  non-convex  or  discontinuous  Pareto  optimal 
fronts. 

The  method  simplifies  the  task  of  the  designer  or  design  team  by  integrating  several  components  so  that 
they  can  focus  on  the  engineering  problem  itself.  The  methods  are  developed  in  a  sequence  of  steps 
consisting  of:  defining  the  requirements,  formulating  the  method,  identifying  several  promising  robust 
analysis  and  optimisation  tools,  the  creation  of  algorithms  and  practical  real-world  problems  in  aeronautics. 

Section  2  of  the  report  describes  the  concept  of  multidisciplinary  design  optimisation;  section  3  the  main 
components  of  the  MDO  methodology.  Section  4  details  the  Evolutionary  Optimiser;  Section  5  details  and 
tests  the  aero-structural  solver  for  two  baseline  designs.  Section  6  details  the  aero-structural  optimisation 
process  and  section  7  presents  the  application  of  the  method  for  two  test  cases  related  to  aero-structural 
UAV  wing  design  optimisation. 
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2  Multidisciplinary  Design  Optimisation 


In  aerospace  engineering  the  designer  is  usually  presented  with  a  problem  which  involves  considering  not 
only  one  objective  but  numerous  objectives  and  multi-physics.  A  systematic  approach  that  accounts  for  the 
coupling  between  the  disciplines  and  variables  regarded  as  Multidisciplinary  Design  Optimisation  (MDO), 
is  required.  In  aircraft  design,  for  example,  the  multi-physics  include  aerodynamics,  structures,  propulsion 
and  control.  These  multi-physics  are  interrelated  and  interdisciplinary  constraints  must  be  satisfied  to  solve 
the  problem.  The  area  of  MDO  has  matured  itself  as  a  separate  discipline  as  evidenced  in  different  journals, 
specialised  technical  committees,  conferences  and  publications  devoted  solely  to  the  topic. 

The  need  for,  and  benefits  of  MDO  are  clear  given  that  even  a  small  improvement  in  the  performance  of 
the  aircraft  can  have  significant  impact  on  the  overall  operation  of  the  vehicle.  An  MDO  process  also 
allows  an  evaluation  of  the  constraints  on  multiple  disciplines  from  the  early  stages  of  the  design,  thus  the 
expense  of  re-designing  an  aircraft  system  is  reduced.  A  comprehensive  survey  of  MDO  methods,  their 
development  and  limitations  is  provided  by  Sobieski  (Sobieski  and  Haftka  1996).  The  research  has 
classified  the  different  methods  and  highlighted  some  important  needs,  including  a  multi-platform 
operation,  the  use  of  parallel  computation  to  reduce  computational  expense  and  space  visualisation  as  the 
designer  might  be  interested  in  the  space  around  the  optimum  rather  than  the  optimum  itself.  In  most  of  the 
methods  described  in  the  survey,  the  optimisation  algorithms  for  MDO  use  traditional  gradient  methods  for 
the  solution,  but  as  has  been  discussed,  these  methods  have  some  limitations.  Sobieski  and  Haftka 
(Sobieski  and  Haftka  1996)  have  evaluated  recent  developments  in  multidisciplinary  aerospace  design  and 
optimisation  and  identified  several  categories  of  problem  formulations  and  also  two  main  challenges  for 
MDO  -  the  computational  expense  and  organisational  (architectural)  complexity.  A  more  recent  survey  by 
Giesing  and  Barthel  (Giesing  1998)  identifies  several  industrial  applications  and  summarises  some  of  the 
needs  for  MDO.  Most  of  the  applications  are  related  to  detailed  design;  few  applications  are  developed  for 
conceptual  or  preliminary  MDO  studies.  On  the  classification  of  needs,  their  research  also  describes  how  a 
MDO  framework  should  be  flexible  to  accept  whatever  functions  is  needed  and  should  address  the  issue  of 
low  -  and  high  -  fidelity  models,  but  not  to  compromise  the  optimisation.  It  also  points  out  the  need  for 
efficient  models  that  describe  the  physics  to  keep  computing  time  at  a  reasonable  level.  A  critical  issue 
mentioned  in  the  paper  is  the  need  for  accurate  and  realistic  design  by  identification  of  the  constraints, 
mechanisms  and  underlying  physics  of  the  various  disciplines  involved. 

An  open  issue  with  MDO  studies  is  the  fact  that  many  high-fidelity  processes  such  as  CFD  or  FEA  are 
complex  to  couple  within  MDO  as  many  of  them  are  not  automated,  robust  or  fast  enough.  Also,  the  need 
for  approximation  techniques,  such  as  Response  Surface  Methodology  (RSM)  (Giunta  1997)  or  Design  and 
Analysis  of  Computer  Experiments  (DACE)  (Sacks,  Welch  et  al.  1989),  arises  because  of  the 
computational  expense  of  using  an  analysis  code  for  all  the  evaluations  during  the  optimisation  process  and 
also  from  the  fact  that  some  analysis  codes  cannot  be  directly  integrated  with  the  MDO  architecture. 

In  MDO  framework  architectures,  there  is  a  requirement  for  more  efficient,  robust  flexible  framework 
architectures  and  methods  with  industrial  codes.  These  codes  should  be  easily  coupled  and  reconfigurable, 
and  adaptable  to  commercial  solvers.  One  of  the  problems  with  current  MDO  architectures  is  that  these  are 
usually  developed  at  universities  and  in  industries  with  restrictions  on  their  use,  and  sometimes  these  are 
specialised  and  difficult  to  generalise  for  other  applications. 

When  implementing  the  optimisation  tools  for  MDO,  it  is  important  to  determine  the  presence  of  multi¬ 
modalities,  non-linearities  and  multiple  objectives  that  might  cause  a  traditional  deterministic  method  to 
fail.  The  computational  cost  of  calculating  the  gradients  and  the  presence  of  multiple  disciplines  is  also  an 
important  consideration. 

The  continuing  challenge  has  been  on  developing  and  improving  numerical  optimisation  techniques  and 
enhancing  their  speed  and  robustness  for  their  use  within  MDO.  One  of  the  emerging  optimisation 
techniques  for  MDO  is  Evolutionary  Algorithms  (EAs),  but  they  have  found  limited  application  in  MDO 
due  to  the  computational  burden  associated  with  them. 
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Giesing  and  Barthel  (Giesing  1998)  also  presented  a  short  discussion  on  supporting  design  space  search 
methods  such  as  EAs,  and  explained  how  they  are  gaining  popularity  for  MDO  because  they  are  simple  to 
couple  with  analysis  modules  and  do  not  incur  the  cost  of  computing  derivatives. 

One  example  of  MDO  is  the  minimisation  of  wing  drag  and  weight  in  a  two-discipline  problem:  namely 
aerodynamics  and  structures.  The  design  variables  can  be  wing  sweep,  wing  span,  twist  or  taper.  With  an 
integrated  approach  the  optimiser  sends  the  values  of  these  variables  to  an  integral  system  of  equations  that 
represent  the  aerodynamic  and  structural  analysis  codes.  These  calculations  are  iteratively  performed  to 
conform  and  converge  on  each  discipline  to  a  consistent  solution.  The  objective  function  and  constraints 
are  then  evaluated  and  manipulated  by  the  optimiser  to  improve  the  design. 

This  approach  is  conceptually  very  simple.  Once  all  disciplines  are  coupled  to  form  a  single 
multidisciplinary  analysis  module,  the  same  techniques  that  are  used  for  a  single-discipline  optimisation 
may  be  used.  The  disadvantage  of  this  approach  is  that  the  solution  of  a  single  system  could  be  very 
expensive  and  does  not  enable  the  potential  decoupling  of  the  individual  disciplines  into  analysis  modules 
that  can  be  computed  in  parallel. 

As  will  be  illustrated  in  Section  7.1,  the  benefits  of  parallelisation  can  be  exploited  by  using  an  EA  which 
sends  candidate  individuals  to  different  processors  where  they  are  evaluated  by  the  complete  MDO  analysis 
and  returned  to  the  optimiser. 

2.1  Limitations  of  Current  Optimisation  Techniques  for  MDO 

As  discussed  in  the  previous  sections,  most  of  the  optimisation  methods  used  for  MDO  use  traditional 
deterministic  techniques.  While  the  use  of  these  methods  is  largely  successful  and  efficient  in  finding 
optimal  global  solutions,  problems  do  arise.  Gradient-based  methods  usually  work  best  with  unimodal 
functions,  but  their  effectiveness  decreases  with  the  presence  of  local  optima  or  ridges  in  the  fitness 
landscape.  Also,  the  presence  of  numerical  noise  inhibits  the  application  of  many  gradient-based 
optimisation  techniques.  The  numerical  noise  causes  an  inaccurate  calculation  of  gradients,  which  in  turn 
slows  or  prevents  convergence  during  optimisation  as  shown  by  Giunta  (Giunta  1997). 

In  aircraft  design,  for  example,  the  problem  of  numerical  noise  is  of  special  concern  when  an  accurate 
solution  is  sought  through  a  high-fidelity  analysis  but  the  computation  of  gradients  is  complex  and  a  single 
aerodynamic  or  structural  analysis  might  take  several  CPU  hours  on  a  supercomputer. 

An  emerging  optimisation  technique  for  MDO  is  Evolutionary  Algorithms  (EAs);  these  techniques  are 
robust  in  finding  optimal  solutions  for  single-  and  multi-objective  problems,  but  have  found  limited 
applications  to  MDO  due  to  the  computational  burden  associated  with  them.  The  challenge  is  then  to  study, 
develop,  apply  and  improve  the  speed  and  robustness  of  these  methods  so  that  confident  applications  and 
within  an  MDO  operational  framework  is  possible. 

It  is  important  to  highlight  that  in  this  work  the  use  of  EAs  will  be  restricted  to  conceptual  and  preliminary 
MDO  studies  where  the  number  of  variables  is  still  relatively  small,  less  than  a  hundred,  and  where  the  use 
of  EAs  is  still  of  potential  benefit.  On  a  larger  scale,  the  use  of  EAs  can  be  extended  for  an  increased 
number  of  variables  and  coupled  with  other  techniques  such  as  Design  of  Experiments  (DOE). 
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3  Formulation  Of  The  MDO  Problem 

The  general  MDO  is  formulated  to  seek  a  vector  X  of  design  variables  that  minimises  an  objective  function 
F(X)  subject  to  disciplinary  interdisciplinary  and  side  constraints  according  to  (Xj)L  <X<  (X)Ut  i  = 
where  N  is  the  number  of  design  variables  and  (Xj,  and  (X)v  are  the  lower  and  upper  bounds  placed  upon 
the  design  variable  Xt.  Two  different  classes  of  design  variables  are  specified: 

•  Aerodynamics 

•  Structures 


Aero  Structural  Optimiser 


Figure  2:  Overall  optimisation  process 


3.1  Method 

The  method  couples  the  proven  Evolutionary  Optimiser  (HAPMOEA)  and  an  Aero- Structural  solver.  The 
Aero- Structural  solver  integrates  two  commercial  high  fidelity  analysis  tools  for  FEA  and  CFD  analysis 
using  MATLAB  as  the  controlling  operating  system. 
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The  optimisation  process  consists  of  the  eight  main  steps  as  illustrated  in  Figure  2.  In  the  first  step,  the 
design  variables,  design  constraints,  flow  conditions  and  fitness  functions  are  defined.  The  second  step 
consists  of  generating  an  initial  population  of  wing  geometries  at  random.  Then  while  the  stopping 
condition  has  not  been  reached,  the  optimiser  generates  new  candidate  geometries  in  steps  three  and  four 
respectively.  In  step  five,  each  candidate  geometry  is  evaluated  by  the  aero-structural  solver.  These 
analyses  provide  the  necessary  information  to  compute  the  fitness  function  or  functions  in  step  six.  If  the 
problem  is  multi-objective  the  optimiser  computes  the  Pareto  fronts  in  step  seven.  The  optimisation 
terminates  if  the  stopping  condition  has  been  met  in  step  eight. 
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4  Evolutionary  Optimiser 
4. 1  EA  Fundamentals 

Evolutionary  algorithms  are  design  and  optimisation  algorithms  which  mimic  the  natural  process  of 
'survival  of  the  fittest'  within  a  generation  of  candidate  solutions.  Broadly  speaking  EAs  operate  simply 
through  the  iterated  mapping  of  one  population  of  solutions  to  another  population  of  solutions.  This  is  in 
contrast  with  the  conventional  deterministic  search  techniques  such  as  the  simplex  method,  conjugate 
gradient  method  and  others,  which  proceed  from  one  given  sub-optimal  solution  to  another,  until  a 
converged  optimum  solution  is  reached.  EAs  are  not  deterministic,  so  that  for  identical  problems  and 
identical  starting  conditions,  the  evolution  of  the  solution  will  not  follow  the  same  path  on  repeated 
applications.  It  is  for  this  reason  that  EAs  fall  into  the  category  of  stochastic  (randomised)  optimisation 
methods.  Some  other  stochastic  methods  that  are  used  are  the  Monte-Carlo  approach,  the  directed  random 
walk  and  simulated  annealing.  The  process  of  evolution  in  EAs  is  of  course  not  completely  random  because 
in  this  case  the  performance  of  the  algorithm  would  be  no  better  than  simple  guessing,  and  at  worst  would 
be  equivalent  to  complete  enumeration  of  the  parameter  search  space.  Evolutionary  algorithms  work  by 
exploiting  population  statistics  to  some  greater  or  lesser  extent,  so  that  when  newer  individual  solutions  or 
offspring  are  generated  from  parents,  some  will  have  inferior  characteristics  and  some  will  have  superior 
characteristics.  The  general  working  principles  of  the  iterated  mapping  then  reduces  to  generating  an 
offspring  population,  removing  a  certain  number  of  inferior  individuals,  and  obtaining  the  subsequent 
population.  This  can  be  summarised  as  the  repeated  application  of  two  operators  on  the  population,  the  first 
being  the  variation  operator  (the  generation  of  offspring)  and  the  second,  the  selection  operator  (the 
survival  of  the  fittest)  (Goldberg  1 989) 

The  origin  of  evolutionary  algorithms  for  parameter  optimisation  seems  to  have  appeared  independently  in 
two  separate  streams,  the  Genetic  Algorithm  (GA)  and  the  Evolution  Strategy  (ES).  The  Genetic  Algorithm 
was  founded  on  principles  developed  by  Holland  (Holland  1975),  and  has  been  extensively  researched  and 
applications  developed.  It  is  generally  accepted  however,  that  the  modern  GA  was  placed  on  its  strong 
foundation  in  optimisation  research  by  Goldberg  (Goldberg  1989).  His  initial  applications  of  the  GA  were 
in  real-world  topics  such  as  gas  pipeline  control.  The  original  GA  technique  revolved  around  a  single 
binary  string  (or  base  -  2)  encoding  of  the  chromosomes,  which  is  the  genetic  material  each  individual 
carries.  The  binary  coded  GA's  variation  operator  is  comprised  of  two  parts,  crossover  and  mutation. 
Crossover  interchanges  portions  of  parental  chromosomes  while  mutation  involves  the  random  switching  of 
letters  in  the  chromosome.  The  selection  operator  has  taken  many  forms,  the  most  basic  being  the 
stochastic  fitness-proportionate  (or  roulette  wheel)  method  (Goldberg  1989). 

Evolutionary  techniques  other  than  GA’s  and  ES’s  exist,  such  as  Evolutionary  Programming  (EP)  and 
Genetic  Programming  (GP).  Evolutionary  Programming  has  been  applied  to  real  coded  optimisation 
problems  (Koza  1992)  but  has  not  seen  widespread  use  in  this  field.  For  comparisons  between  ES  and  EP 
refer  to  Back  (Back,  Rudolph  et  al.  1993).  Genetic  Programming  (Koza  1992)  is  devoted  to  the  generation 
of  computer  programs  rather  than  number  sets  as  the  solution  to  a  given  problem.  Together  with  GA’s  and 
ES’s  these  methodologies  form  the  basic  four  'schools'  of  evolutionary  algorithms  (GA,  ES,  EP  and  GP). 

Evolution  Strategies  were  first  developed  at  the  Technical  University  of  Berlin,  (Back,  Rudolph  et  al. 
1993).  The  first  algorithm  worked  using  only  two  individuals,  one  parent  and  one  offspring.  Each 
individual  was  real  coded;  each  problem  variable  was  assigned  a  floating  point  value  in  the  chromosome. 
The  variation  operator  involved  applying  a  random  mutation  to  each  floating  point  value  in  the  parental 
chromosome  to  arrive  at  the  offspring  individual.  The  selection  operator  was  entirely  deterministic,  and 
was  simply  the  result  of  a  competition  between  parent  and  offspring  to  determine  which  one  remained.  In 
the  standard  nomenclature  this  strategy  is  denoted  the  (1+1)  ES,  the  first  digit  indicating  the  number  of 
parents,  the  '  +  '  indicating  competition  between  parents  and  offspring  and  the  final  digit  indicating  the 
number  of  offspring’s.  From  the  beginning  the  ES  has  been  designed  almost  exclusively  with  real  coding  in 
mind,  as  opposed  to  original  GA  variants  where  real  parameter  optimisation  comes  about  by  the  piecewise 
interpretation  of  the  binary  chromosome  associated  with  each  individual.  An  evolution  strategy  would 
therefore  seem  a  logical  starting  point  for  evolutionary  optimisation  using  real  coded  problem  variables. 
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Subsequent  developments  in  ES’s  introduced  multi-membered  populations,  and  the  first  algorithm  of  this 
type  was  the  (p+1)  ES  (Back  and  Schwefel  1995).  This  worked  by  applying  some  variation  operator  to  p 
parents  to  produce  a  single  offspring.  The  offspring  is  selected  by  determining  whether  it  is  better  than  the 
worst  member  of  p~,  and  if  so  it  replaces  the  worst  member  p.  Both  the  (1+1)  ES  and  the  (p+1)  ES  used 
deterministic  control  of  the  mutation  size  (variations  applied  to  design  variables)  which  were  normally 
distributed  when  applied  to  real  coded  problems. 

A  pseudo  code  of  a  canonical  evolution  strategy  is  illustrated  in  Figure  3.  A  population  (ju0)  is  initialised 
and  then  evaluated.  For  a  number  of  generations  (g),  and  while  a  stopping  condition  (maximum  number  of 
function  evaluation  or  target  fitness  value)  is  not  yet  met,  offspring  (l8+i)  go  recursively  through  the 
process  of  recombination,  mutation,  evaluation  and  selection. 


Initialise : 

init{jUo ) 

Evaluate : 

/(a) 

9=0 

while  stopping  condition  not  met. 

Recombine:  =  rCCoifU S  ) 

Mutate : 

X^=mut{xfl) 

Evaluate : 

^+I=/fe+1) 

Select : 

jLl8+l  =  Sel^jLl  +  X)  (Plus  strategy)  or, 

(l **  =  Sel(X)  (Comma  strategy) 

9=9+1 

loop 

Figure  3:  Label  Canonical  Evolution  Strategy 


The  recent  developments  in  both  GAs  and  ESs  have  greatly  modified  their  variation  and  selection 
operators,  to  the  point  where  it  is  not  clear  whether  such  a  nomenclature  division  is  nowadays  particularly 
justified.  The  main  difference  that  exists  between  them  today  is  still  the  predominance  of  adaptive 
mutations  in  ES’s,  which  have  made  them  very  attractive  for  real  coded  optimisation,  although  GA  research 
has  produced  some  related  concepts  (Collard  and  Escazut  1995;  Sefioui,  Periaux  et  al.  1996;  Thierens 
2002) 

Some  of  the  advantages  of  EAs  are  that  they  require  no  derivatives  or  gradients  of  the  objective  function, 
have  the  capability  of  finding  globally  optimum  solutions  amongst  many  local  optima,  are  easily  executed 
in  parallel  and  can  be  adapted  to  arbitrary  solver  codes  without  major  modification.  Another  major 
advantage  of  EAs  is  that  they  can  tackle  multi-objective  problems  directly.  EAs  have  been  successfully 
applied  to  different  aeronautical  design  problems  and  there  have  been  different  efforts  to  explore  the 
benefits  and  capabilities  of  EAs  for  MDO,  aircraft,  wing,  aerofoil  and  rotor  blade  design  (Makinen  1998; 
Obayashi  1998;  Sasaki  2004;  Gonzalez  2005). 

4.2  The  Development  of  Evolutionary  Algorithms  for  Design  and 
Optimisation  in  Aeronautics 

The  potential  benefits  of  EAs  for  optimisation  problems  in  engineering  have  been  recognised  for  some 
time.  Different  studies  explore  the  potential  benefits  of  EAs  for  wing  design  and  optimisation  (Obayashi 
1998),  (Oyama,  Liou  et  al.  September,  2002),  (Oyama  2000).  Obayashi,  for  example,  has  applied  EAs  for 
several  wing  platform  design  problems 
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In  these  studies,  different  niching  and  elitist  models  are  applied  to  a  Multi  Objective  Genetic  Algorithm 
(MOGA)  to  find  optimal  Pareto  fronts  for  transonic  and  supersonic  wing  design.  The  transonic  case 
considered  three  (3)  objectives:  minimise  aerodynamic  drag,  minimise  aircraft  weight  and  maximise  fuel 
weight  stored  in  the  wing.  The  constraints  imposed  were  lift  greater  than  given  aircraft  weight  and 
structural  strength  greater  than  aerodynamic  loads. 

Other  applications  of  EAs  for  wing  design  include  the  work  by  Obayashi  (Obayashi  1998)  ,  (Anderson 
1996;  Takahashi  1997),  (Gonzalez,  Whitney  et  al.  2004).  Gonzalez,  Whitney  et  al.  (Gonzalez,  Whitney  et 
al.  2004  )  for  example  studied  and  illustrated  the  application  of  a  hierarchical  topology  of  EAs  for  multi¬ 
objective  multidisciplinary  wing  design  using  low  fidelity  analysis  tools.  The  main  results  of  these  studies 
indicate  the  importance  of  variable  fidelity  models,  the  broad  applicability  and  the  ability  of  EAs  to  find 
optimal  Pareto  solutions  for  three-dimensional  applications  and  problems  with  more  than  two  objectives. 

Crispin  (Crispin  1994  )  applied  GAs  for  aircraft  conceptual  and  preliminary  design  and  found  it  useful  to 
obtain  reasonable  and  feasible  solutions.  Crossley  et  al  (Crossley  1996)  applied  a  GA  for  helicopter 
conceptual  design  and  were  able  to  show  the  effectiveness  of  GAs  and  obtain  optimal  feasible 
configurations.  One  of  the  results  of  his  work  was  that  the  use  of  parametric  variations  conducted  by  GAs 
can  significantly  reduce  the  amount  of  time  and  money  in  the  early  stages  of  aircraft  design. 

Ruben  et  al.  (Ruben  and  .  2000)  conducted  some  work  on  GAs  for  aircraft  conceptual  design  and  found  a 
5%  weight  savings  when  compared  to  a  conventional  design.  Ali  and  Behdinan  (Ali  2002)  applied  GAs  to 
determine  an  optimal  combination  of  design  variables  for  a  medium-size  transport  aircraft.  They  studied 
different  selection  and  crossover  strategies  and  indicated  that  with  a  GA  approach  it  was  possible  to 
generate  feasible  and  efficient  conceptual  designs.  In  these  studies,  the  authors  also  highlighted  the 
effectiveness  and  importance  of  EAs  in  saving  money  in  the  initial  stages  of  the  design  process. 

Parmee  and  Watson  (Parmee  and  Watson  1999)  proposed  a  preliminary  airframe  design  using  co¬ 
evolutionary  multi-objective  genetic  algorithms.  Their  algorithm  was  able  to  find  local  objective  optimal 
solutions  after  a  few  generations  and  identify  paths  to  trace  the  trade-off  surface  to  some  extent.  This 
research  also  mentions  that  on-line  sensitivity  analysis  has  a  role  to  play  as  the  number  of  objectives 
increases  and  suggests  that  quicker,  less  detailed  runs  can  easily  be  achieved  using  smaller  population  sizes. 

Other  applications  of  EAs  for  aircraft  design  include  the  work  by  Cvetkovick  et  al.  (Cvetkovic  2000)  and 
Raymer  (Raymer  2002;  Raymer  2003).  Traditional  transport  or  commercial  aircraft  have  been  studied  and 
optimised  using  EAs.  A  few  studies  on  the  benefits  of  EAs  in  exploring  large  variations  in  the  design  space 
have  been  conducted  for  novel,  non-notional  configurations  such  as  Unmanned  Aerial  Vehicles 
(UAVs/UCAVs)  and  Micro  Aircraft  Vehicles  (Ng  2002),  (Gonzalez,  Whitney  et  al.  2004),  (Gonzalez, 
Whitney  et  al.  2004  ),  (Whitney  2003). 

All  these  works  provide  an  indication  of  the  broad  applicability  and  robustness  of  Evolutionary  Algorithms 
to  find  optimal  solutions  and  how  industry  is  increasingly  applying  them  to  solve  complex  problems  in 
engineering. 

Little  work  has  been  conducted  on  comparing  the  application  of  EAs  and  other  methods  for  MDO.  Raymer 
and  Crossley  (Raymer  2002;  Raymer  2003)  applied  and  compared  different  MDO  methods,  Monte  Carlo, 
Random  Walk,  Simulated  Annealing,  GAs  and  orthogonal  steepest  descent  search  to  enhance  aircraft 
conceptual  design  and  MDO.  In  his  research  Raymer  applied  these  methods  to  four  aircraft  concepts:  a 
fighter,  a  commercial  airliner,  an  asymmetrical  light  twin  and  a  tactical  UAV  showing  a  broad  applicability 
of  GAs.  One  of  the  conclusions  of  his  work  is  that  aircraft  conceptual  design  can  be  improved  by  proper 
application  of  optimisation  methods  for  MDO.  He  found  that  the  proper  selection  of  a  technique  can  reduce 
the  weight  and  cost  of  an  aircraft  concept  by  minor  changes  in  the  design  variables.  His  results  also  indicate 
that  the  orthogonal  steepest  descent  method  provides  a  slightly  better  result  with  the  same  number  of 
function  evaluations,  but,  as  the  number  of  variables  is  increased,  the  evolutionary  methods  seem  to  be 
superior.  It  should  be  noted  that  Raymer  limited  his  research  to  single-objective  problems,  used  only  one 
type  of  fidelity  model  for  the  aircraft  analysis  and  limited  his  research  to  seven  design  variables  and  the 
inclusion  of  any  propulsion  variables  other  than  engine  size  via  Thmst/Weight. 
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4.3  Advantages  and  Limitations  of  Traditional  EAs  for 

Aeronautical  Problems 

Even  though  there  are  definite  advantages  to  using  EAs  they  have  not  seen  a  widespread  use  in  industry  or 
for  multidisciplinary  design  optimisation  applications.  The  main  reason  for  this  is  the  computational 
expense  involved,  as  they  require  more  function  evaluations  than  gradient-based  methods.  Therefore,  the 
continuing  challenge  has  been  to  improve  their  performance  and  develop  new  numerical  techniques. 

It  is  clear  that  a  possible  and  realistic  avenue  is  the  application  of  robust  and  efficient  evolutionary  methods 
for  MDO.  One  of  the  viable  alternatives  is  the  Hierarchical  Asynchronous  Parallel  Multi-Objective 
Evolutionary  Algorithms  (HAPMOEA)  described  in  the  next  section.  These  algorithms  have  proven  to  be 
robust  and  will  be  extended  in  this  research  for  multi-objective  and  multidisciplinary  design  optimisation 
problems  in  aeronautics. 

4.4  Hierarchical  Asynchronous  Parallel  Multi-objective 

Evolutionary  Algorithms 

One  drawback  of  EAs  is  that  they  are  slow  to  converge  in  that  they  require  a  large  number  of  function 
evaluations  to  refine  the  global  optimum  position  once  found.  Gradient  based  optimisation  techniques  do 
not  suffer  from  this  problem.  Furthermore,  a  EAs  performance  decreases  as  the  number  of  variables 
increase.  Hence  the  continuing  challenge  with  EAs  has  been  to  develop  robust  and  fast  numerical 
techniques  to  overcome  these  known  difficulties  and  facilitate  the  complex  task  of  design  and  optimisation 
in  aeronautics.  In  this  direction  we  developed  and  use  the  Hierarchical  Asynchronous  Parallel  Evolutionary 
Algorithm  (HAPMOEA)  approach  (Whitney,  Sefrioui  et  al.  2002;  Gonzalez  2005;  Gonzalez  2005)  with 
some  extensions  for  Multi-Disciplinary  and  Multi-Objective  analysis  introduced  recently.  The  foundations 
of  the  algorithm  lie  upon  traditional  evolution  strategies  and  incorporate  the  concepts  of  covariance  matrix 
adaptation  (CMA)  (Hansen  and  Ostermeier  2001),  multi-objective  optimisation,  hierarchical  topology, 
parallel  computing  (parallel  EAs)  and  asynchronous  evaluation  of  candidate  solutions,  Pareto  tournament 
selection  for  single  and  multi-objective  problems  and  a  constraint  handling  mechanism. 

4.4.1  Multi-Criteria  Optimisation 

Often  aeronautical  design  problems  require  a  simultaneous  optimisation  of  conflicting  objectives  and 
associated  number  of  constraints.  Unlike  single  objective  optimisation  problems,  the  solution  is  a  set  of 
points  known  as  Pareto  optimal  set.  Solutions  are  compared  to  other  solutions  using  the  concept  of  Pareto 
dominance.  A  multi-criteria  optimisation  problem  can  be  formulated  as: 

Maximise  /  Minimise 

. ’N 

Subject  to  constraints : 

g/(*)  =  °  7=1, M 

hk(x)<0  k- 1, ,M 

Where  f  are  the  objective  functions,  N  is  the  number  of  objectives;  x  is  an  n  -  dimensional  vector  where  its 
arguments  are  the  decision  variables.  For  a  minimisation  problem,  a  vector  xt  is  said  partially  less  than 
vector  x2  if: 


V,  (v)  A/«:  (a  )  and  (v)<./:  (x2  ) 


In  this  case  the  solution  xj  dominates  the  solution  x2. 
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Figure  4  shows  Pareto  optimality  for  a  two  objective  problem.  The  objective  of  the  optimisation  is  then  to 
provide  a  set  of  Pareto  optimal  solutions  that  represent  a  trade-off  of  information  amongst  the  objectives. 


Figure  4:  Pareto  Optimality 


As  EAs  evaluate  multiple  populations  of  solution  points,  they  are  capable  of  finding  a  number  of  solutions 
in  a  Pareto  set.  Pareto  selection  ranks  the  population  and  selects  the  non-dominated  individuals  for  the 
Pareto  fronts.  An  Evolutionary  Algorithm  that  has  capabilities  for  Multi-Objective  optimisation  is  termed 
Multi-Objective  Evolution  Algorithms  (MOEAs).  Theory  and  applications  of  MOEAs  can  be  found  in  Deb 
(Deb  2003)  and  Coello-Coello  et  al  (Coello,  Veldhuizen  et  al.  2002). 

4.4.2  Hierarchical  Population  Topology 

A  hierarchical  population  topology,  when  integrated  into  an  evolution  algorithm,  means  that  a  number  of 
separate  populations  are  established  in  a  hierarchical  layout  to  solve  the  given  problem,  rather  than  a  single 
‘cure-all1  type  single  population  layout.  This  method  was  proposed  by  Sefrioui  (Sefrioui  and  Periaux  2000) 
and  is  shown  in  Figure  5.  The  bottom  layer  can  be  entirely  devoted  to  exploration,  the  intermediate  layer  is 
a  compromise  between  exploitation  and  exploration  and  the  top  layer  concentrates  on  promising  solutions. 
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4.4.3  Parallelisation  Of  The  Algorithm 

The  algorithm  is  similar  to  hybrid  parallel  Multi-Objective  Evolutionary  Algorithms  (pMOEA)  (Cantu-Paz 
2000;  Veldhuizen,  Zydallis  et  al.  2003);  it  uses  a  master  slave  approach  but  incorporates  the  concept  of 
isolation  and  migration  through  a  hierarchical  binary  tree  starcture  where  each  level  has  different 
parameters  for  the  MOEAs  .  The  parallel  environment  used  in  this  work  is  a  cluster  of  PCs,  wherein  the 
master  carries  on  the  optimisation  process  while  remote  nodes  perform  the  analysis  (CFD,  FEA).  The 
message-passing  model  used  is  the  Parallel  Virtual  Machine  (PVM)  (Geist,  Beguelin  et  al.  1994).  The 
algorithm  has  been  tested  in  a  cluster  of  heterogeneous  CPUs  available  at  the  School  of  Aerospace, 
Mechanical  and  Mechatronic  Engineering  at  the  University  of  Sydney  and  with  different  RAMs 
frequencies,  caches,  memory  access  times,  storage  capabilities  and  communication  attributes.  The  cluster 
can  be  configured  with  up  to  18  machines  with  performances  varying  between  2.0  and  3.4  GHz. 

4.4.4  Asynchronous  Evaluation 

When  considering  the  solution  to  Multi-objective  and  Multidisciplinary  Optimisation  cases,  several 
problems  arise  as  many  methods  of  solution  used  in  engineering  today  may  take  different  times  to  complete 
their  operation.  The  classic  example  of  this  is  the  modern  CFD  solver.  With  a  typical  industrial  code  used 
for  external  aerodynamic  analysis  of  airplanes,  the  time  for  the  residual  of  the  solution  to  converge  to  a 
specified  level  (either  machine  zero  or  an  arbitrarily  selected  higher  value)  can  vary  over  a  significant  range 
depending  on  the  geometry,  boundary  condition,  the  algorithm  and  CPU  used.  The  previous  generation  of 
evolutionary  algorithms  have  mostly  used  a  generation-based  approach  and  so  are  the  traditional  genetic 
algorithm  and  evolution  strategy.  A  difficulty  with  generational  models  is  that  they  may  create  an 
unnecessary  bottleneck  when  used  on  parallel  computers.  If  the  population  size  is  approximately  equal  to 
the  number  of  processors,  and  most  of  the  candidate  offspring’s  that  are  sent  for  solution  can  be 
successfully  evaluated,  it  may  so  happen  that  some  processors  will  complete  their  task  quicker  with  the 
remainder  taking  more  time.  With  a  generational  approach,  those  processors  that  have  already  completed 
their  solutions  will  remain  idle  until  all  processors  have  completed  their  work. 
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The  approach  used  here  is  to  ignore  any  concept  of  a  generation  based  solution  and  is  similar  to  work  by 
Wakunda  and  Zell  (Wakunda  and  Zell  2002)  and  other  non-generational  approaches.  The  difference  is  that 
the  selection  operator  couples  a  one-by-one  (steady-state)  function  evaluation  with  a  direct  multi-objective 
fitness  criterion.  Whilst  a  parent  population  exists,  offspring  are  not  sent  as  a  complete  ‘block’  space  to  the 
parallel  slaves  for  solution.  Instead  one  candidate  is  generated  at  a  time,  and  is  sent  to  any  idle  processor 
where  it  is  evaluated  at  its  own  speed.  When  candidates  have  been  evaluated,  they  are  returned  to  the 
optimiser  and  either  accepted  by  insertion  into  the  main  population  or  rejected.  This  requires  a  new 
selection  operator  because  the  offspring  cannot  now  be  compared  one  against  the  other,  or  even  against  the 
main  population  due  to  the  variable-time  evaluation.  The  recently  evaluated  offspring  are  compared  against 
a  previously  established  rolling-benchmark  and  if  successful,  replace  (according  to  some  rule)  a  pre¬ 
existing  individual  in  the  population.  This  benchmarking  is  implemented  via  a  separate  evaluation  buffer 
( B ),  which  provides  a  statistical  ‘background  check'  on  the  comparative  fitness  of  the  solution.  The  length 
of  the  buffer  should  represent  a  reasonable  statistical  sample  size,  but  need  not  be  too  large;  approximately 
twice  the  population  size  is  more  than  ample.  When  an  individual  has  had  a  fitness  assigned,  it  is  then 
compared  to  past  individuals  (both  accepted  and  rejected)  to  determine  whether  or  not  it  should  be  inserted 
into  the  main  population.  If  it  is  to  be  accepted,  then  some  replacement  strategy  is  invoked  and  it  replaces  a 
member  of  the  main  population.  The  replace-worst-always  method  is  exclusively  used  in  this  work.  A 
schematic  of  the  asynchronous  and  parallelization  approach  is  shown  in  Figure  5 


Figure  6:  Parallel  Computing  and 
Asynchronous  Evaluation 


4.4.5  Pareto  Tournament  Selection 

An  on-the-fly  selection  operator  is  implemented  by  means  of  a  Pareto  tournament  selection  operator.  To 
implement  an  optimisation  algorithm  that  is  equally  applicable  to  both  single  and  multi-objective  problems, 
a  suitable  selection  operator  capable  of  handling  either  situation  must  be  developed.  An  extension  of  the 
standard  tournament  operator  popular  in  many  approaches  is  proposed  (Goldberg  1989;  Michalewicz 
1992).  Most  evolutionary  algorithms  configured  for  multi-objective  optimisation  currently  use  the  non- 
dominated  sorting  approach.  This  is  a  straightforward  way  to  adapt  an  algorithm  that  is  designed  as  a  single 
objective  optimiser  into  a  multi-objective  optimiser,  and  is  used  by  many  researchers  (Michalewicz  1992; 
Deb  2003)].  The  problem  with  sorting  approaches  is  that  the  method  is  not  a  fully  integrated  one.  Briefly,  a 
sorting  method  works  by  computing  the  set  of  non-dominated  solutions  amongst  a  large  statistical  sampling 
(either  a  large  population  or  previous  data),  and  assigning  these  solutions  a  rank  one.  Then  ignoring  these 
points,  the  process  is  repeated  until  a  ‘second’  Pareto  front  is  found,  and  this  is  assigned  rank  two.  This 
process  continues  until  all  points  are  ranked,  and  then  the  value  of  the  rank  is  assigned  to  the  individual  as  a 
new  single  objective  fitness.  A  problem  arises  now  on  whether  it  is  fair  to  assign  individuals  in  the  second 
rank  numerically  half  the  fitness  of  the  first,  and  whether  the  third  rank  deserves  a  third  of  the  fitness  of  the 
first.  This  poses  a  dilemma  regarding  the  level  of  equality  present  amongst  the  solutions,  as  often  solutions 
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with  excellent  information  may  lie  adjacent  to,  but  not  in,  rank  one.  To  solve  this  ‘artificial  scaling’ 
problem,  it  is  possible  to  introduce  scaling,  sharing  and  niching  schemes,  however  all  of  these  require 
problem-specific  parameters  or  knowledge,  even  in  adaptive  approaches.  It  is  of  course  always  desirable  to 
compose  an  algorithm  that  does  not  introduce  such  unnecessary  parameters. 

The  current  operator  is  a  novel  approach  in  that  it  requires  no  additional  ‘tuning’  parameters,  works 
seamlessly  with  the  asynchronous  selection  buffer  (B),  and  is  very  easy  to  encode.  As  illustrated  in  Figure  7 
to  determine  whether  a  new  individual  x  is  to  be  accepted  into  the  main  population,  we  compare  it  with  the 
selection  buffer  by  assembling  a  small  subset  of  the  buffer  called  the  tournament  Q  =  [q1,q2,...,q„\  We 
assemble  Q  by  selecting  individuals  from  the  buffer,  exclusively  at  random,  until  it  is  full.  We  then  simply 
ensure  that  the  new  individual  is  not  dominated  by  any  in  the  tournament.  If  this  is  the  case,  then  it  is 
immediately  accepted,  and  is  inserted  according  to  the  replacement  rules.  The  only  parameter  that  needs  to 
be  determined  in  advance  is  the  tournament  size,  a  parameter  that  would  exist  in  a  single  objective 
optimisation  anyway.  Selection  of  this  parameter  requires  a  small  amount  of  problem  specific  knowledge, 
and  should  vary  between  Q  =  B/2  (strong  selective  pressure)  and  Q  =  B/6  (weak  selective  pressure).  The 
optimiser  is  not  overly  sensitive  to  this  value,  provided  the  user  errs  on  the  side  of  weak  selective  pressure 
(smaller  tournaments)  in  the  absence  of  better  information.  The  egalitarian  approach  to  the  tournament  (by 
selecting  individuals  at  random)  ensures  good  diversity  amongst  the  selected  individuals;  no  niching  or 
forced  separation  of  individuals  has  been  found  necessary.  It  can  also  be  seen  that  in  the  event  the  fitness 
vectors  have  only  one  element  (a  single  objective  optimisation),  this  operator  simplifies  to  the  standard 
tournament  selection  operator  (Goldberg  1989;  Michalewicz  1992). 


Population 

Asynchronous  Buffer  (5) 


Tournament  Q 

Evaluate  x  x 
If  ^  not  dominated 


-B<0<-B 

6  2 


Figure  7:  Pareto  Tournament 
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5  Aero-Structural  Analysis 


5.1  Introduction 

This  section  describes  the  aero-structural  solver  in  detail.  Sectio  5.1  describes  the  aero-structural  program 
layout,  section  5.2  details  the  validation  of  the  aerodynamic  solver  while  section  5.3  describes  the  structural 
model  and  analysis.  Section  5.4  details  the  aero-structural  analysis  and  results  for  two  baseline  geometries. 

5.2  Aero-Structural  Program  Layout 

A  program  to  perform  the  aero-structural  optimization  was  written  in  Matlab®.  The  aero  structural  solver 
program  integrates  two  commercial  analysis  tools  for  Finite  Element  Analysis  (FEA)  and  Computational 
Fluid  Dynamics  (CFD)  namely  MSC.Nastran®,  developed  by  MSC  Software,  PANAIR,  a  high  order 
Mornio  class  panel  method  developed  by  Boeing  and,  FRICTION  to  calculate  the  form  and  frictional  drag 
about  a  candidate  wing.  The  entire  aero-structural  program  is  controlled  through  a  Matlab  "  script  file.  This 
allows  for  an  easy  coupling  of  the  different  required  programs  as  one  continuous  Matlab®  data  stmcture  can 
be  utilized  to  define  all  the  information  passed  between  programs.  A  flow  diagram  of  this  process  is 
illustrated  in  Figure  8.  A  description  of  the  aerodynamic  and  stmctural  analysis  tools  is  given,  followed  by 
a  general  description  of  the  validation  cases. 

The  aero-structural  program  works  in  a  very  structured  way  as  shown  in  Figure  8.  On  execution  of  the 
code,  the  workspace  is  cleared  of  any  existing  variables  and  any  output  files  from  previous  program  runs 
which  may  have  existed  in  the  working  directory.  The  program  loads  in  any  user  or  optimiser  specified 
inputs  through  the  calling  of  the  user  generated  input  files. 

Once  loaded,  the  Aero-Structural  program  performs  a  number  of  geometry  checks  on  the  candidate  wing 
and  then  creates  the  wing  shape  from  the  user  specified  inputs.  From  this  output  the  FRICTION,  PANAIR 
and  MSC.Nastran"  input  files  are  created  and  executed.  After  running  the  PANAIR  simulation,  the 
aerodynamic  loads  are  read  in  from  the  PANAIR  output  files  and  mapped  to  the  MSC.Nastran®  simulation. 

Once  again  the  candidate  wing  is  tested  but  this  time  it  is  to  see  whether  or  not  the  wing  violates  any 
penalty  function  to  do  with  the  aerodynamics  or  stmctures.  Thereafter  an  output  file  containing  the 
simulation  results  is  produced. 
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Figure  8:  Aero-Structural  Analysis  Program  Layout 


If  an  unrecoverable  error  is  encountered  anywhere  in  the  solution  sequence,  the  aero-structural  solver  is 
able  to  recover  and  generate  a  simulation  output  file  which  indicates  to  the  evolutionary  optimiser  that  the 
solution  failed. 

5.1.1  Input  Files 

The  Aero-Structural  solver  loads  in  information  describing  the  wing  from  two  user  defined  sources.  The 
first  source  is  through  the  variables  defined  in  the  ‘ UserWing.m ’  and  “User Auxiliary Data.m”  files.  This  is 
the  most  common  and  simplest  method  used  when  only  performing  a  single  run  of  the  Aero- Structural 
solver.  When  coupled  to  the  Evolutionary  Optimiser,  a  different  method  is  used. 

The  loading  in  of  Evolutionary  Optimiser  input  variables  to  PANAIR  is  performed  by  a  number  of 
functions  written  in  Matlab".  The  functions  have  the  ability  to  load  in  the  variables.  Four  output  files  are 
used  depending  on  the  types  of  variables  being  passed  to  the  aero-structural  solver.  Each  file  contains 
different  information  relating  to  the  wing.  The  four  different  types  of  input  files  are  used  by  the  Aero- 
Structural  Solver  are: 
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•  Planform  data  ( planform  Variables. XX) 

•  Auxiliary  data  (auxiliatyVariables.XX) 

•  Aerodynamics  ( a erofo ils Da ta.XX) 

•  Structural  data  ( structuralVariables.XX) 

These  files  are  simple  comma  separated  value  files  describing  each  variable  passed  to  the  optimiser. 

Planform  data  includes  crank  locations,  sweep  angles,  etc.  Auxiliary  data  handles  the  different  flow 
conditions  such  as  Mach  number  and  angle  of  attack.  Aerodynamics  simply  contains  the  different  aerofoil 
shapes  as  x  and  y  columns  corresponding  to  the  aerofoil  sections  at  the  root,  tip  and  crank  positions.  The 
structural  data  has  the  thickness  of  the  different  elements  making  up  the  aircraft  wing. 

5.1.2  Output  Files 

The  Aero- Structural  Solver  creates  two  output  files.  The  first  file  contains  a  simple  indication  as  to  whether 
or  not  the  Aero-Structural  Analysis  Program  simulation  had  terminated  successfully,  the  second,  if  the 
termination  was  successful,  contains  the  aerodynamic  and  structural  results  from  the  simulation. 

To  indicate  whether  or  not  the  simulation  had  terminated  successfully  (there  were  no  critical  errors 
requiring  Matlab®  to  unexpectedly  return  command  to  the  command  prompt),  a  convergence  file 
‘convergence.txt’  is  created  with  either  a  zero  (0),  (no  errors  and  hence  the  results  file  was  created)  or  a  one 
(1)  to  indicate  errors  within  the  program. 

If  the  aero-structural  solver  is  able  to  mn  all  the  required  simulations  on  a  candidate  wing,  a  results  file  is 
written.  This  results  file  has  a  simple  header  which  is  a  legend  into  the  remaining  lines  in  the  file.  The 
legend  describes  what  the  value  in  the  file  corresponds  to,  be  it  the  induced  drag  calculated  by  PANAIR, 
friction  drag  calculated  by  FRICTION,  penalties  the  wing  has  been  penalized  with  or  the  stmctural 
deflections.  If  the  aero-stmctural  is  unable  to  calculate  any  value  due  to  a  recoverable  error  such  as  the 
MSC.Nastran"  simulation  failing,  the  affected  results  are  indicated  by  a  hash  symbol  (#). 

5.1.3  Automatic  Wing  Geometry  Generation 

The  aerodynamic  and  stmctural  wing  models  are  created  though  an  automatic  geometry  modeller 
developed  for  this  task  (‘ MatDefm ’).  The  aerodynamic  model  is  created  first  and  provided  the  body  from 
which  the  stmctural  model  is  created. 

After  all  the  input  data  checks  are  performed,  the  aerodynamic  model  is  built  in  sections.  Each  section  is 
made  up  of  the  area  between  cranks  as  is  shown  in  Figure  9.  The  first  section  constmcted  is  between  the 
root  and  intermediate  crank  and  the  second  section  is  constmcted  between  the  intermediate  and  tip  crank. 

Matrix  operations  are  used  to  define  the  initial  layout  of  the  section  and  take  into  account  sweep  and  taper 
ratios.  All  sweep  angles  are  defined  relative  to  the  global  XY  frame  defined  according  to  normal  aircraft 
convention.  The  aerofoil  shapes  at  the  intermediate  rib  locations  are  determined  through  a  linear 
interpolation  of  the  aerofoils  defined  at  each  section  boundary. 
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Root  Crank  Intermediate  Crank 


Figure  9:  Plan  view  of  a  Simple  two  section  wing  showing  Construction 


Once  all  the  sections  have  been  created,  the  entire  wing  is  defined  through  the  amalgamation  of  all  the 
computed  sections. 

As  mentioned  above,  the  structural  model  is  extracted  from  the  aerodynamic  model.  As  all  the  cranks  are 
defined  according  to  the  internal  structure,  this  operation  is  easily  performed. 

Along  with  the  external  definition,  when  constructing  the  stmctural  model,  the  automatic  geometry 
modeller  computes  the  location  of  all  the  Spar  and  Rib  Caps. 

5.1.4  Automatic  Certification  Of  The  Wing  Geometry 

Sometimes  a  random  combination  of  number  of  cranks,  chord  ratios,  sweep  angles,  etc  or  material 
properties  can  result  in  a  physically  unfeasible  geometry.  The  Aero- Structural  Analysis  Program  therefore 
has  a  number  of  internal  checks  to  validate  the  data  that  is  passed  to  it. 

Structures 

The  material  data  defined  in  the  ‘ MaterialProperties.m ’  file  is  analysed  to  check  if  enough  information  is 
provided.  Furthermore,  the  program  determines  which  variable  is  missing  data  for  certain  crank  locations 
and  either  replicates  data,  or  if  many  variables  are  missing  data,  remove  data  from  variables  such  that  a 
common  denominator  is  found. 

Geometry 

Two  simple  geometry  checks  are  performed  on  the  wing  before  any  aerodynamic  or  structural  analysis  is 
performed.  The  first  check  ascertains  whether  or  not  the  wing  upper  and  lower  surfaces  intersect  each  other 
(Figure  10  a-Crossing),  and  the  second  check  ascertains  if  the  leading  and  trailing  edges  intersect  (Figure 
10  b-S  witching).  If  either  of  these  checks  fails,  the  simulation  is  aborted. 


23 


Once  both  the  data  checks  and  basic  geometry  checks  have  been  performed,  The  Aero-Structural  Analysis 
Program  checks  to  see  if  any  of  the  data  defining  the  internal  layout  of  the  wing  will  cause  errors  later  in 
the  simulation  steps.  Errors  can  be  caused  by  a  wing  rib  being  referenced  by  two  crank  locations.  This  can 
be  caused  by  the  method  in  which  the  ribs  are  linearly  spaced  throughout  the  wing  and  the  crank  locations 
are  rounded  up  or  down  to  the  nearest  rib  as  in  Figure  1 1 .  This  above  mentioned  case  can  cause  modelling 
problems  and  hence  one  of  the  sets  of  data  defining  one  of  the  crank  locations  would  have  to  be  removed. 


This  removal  process  is  carried  out  by  deleting  the  data  defining  the  more  outboard  of  the  two  or  more  data 
sets  and  retaining  the  inboard. 

5.2  Aerodynamic  Analysis 

A  proper  selection  and  validation  of  analysis  tools  is  required  before  coupling  it  with  an  optimisation 
process.  A  flow  solver  should  meet  some  essential  requirements  such  as:  result  accuracy,  computational 
expense  and  robustness. 

It  is  always  desirable  to  use  a  high  fidelity  solver  that  can  account  for  the  flow  complexities.  The  problem 
of  using  a  full  Navier-Stokes  solver  is  the  computational  expense  of  the  solution  as  one  computation  on  a 
full  3  dimensional  wing  might  take  several  hours  on  a  supercomputer.  Therefore,  to  reduce  this  expense  it 
was  necessary  to  find  another  high  fidelity  method  which  could  yield  accurate  results  at  a  minimum 
computational  cost.  It  was  for  the  above  reason  why  the  authors  chose  to  use  a  high-order  Morino  class 
panel  method  code  called  PANAIR  to  solve  for  the  flow  properties  about  the  arbitrary  vehicle  wing. 
PANAIR  only  calculates  the  Induced  drag  component  produced  by  the  wing  as  it  travels  through  the  air. 
Hence  to  correctly  predict  the  total  drag  produced  by  the  wing,  both  the  Form  and  Friction  components 
need  to  be  calculated;  for  this  task  FRICTION  (Mason  2001)  was  used. 

5.2.1  High  Order  Panel  Method  (PANAIR) 

The  higher  order  panel  method  PANAIR  was  developed  by  Boeing  as  A502  to  predict  subsonic  and 
supersonic  linearised  potential  flows  about  an  arbitrary  geometry  (Saaris  1992).  PANAIR  differs  from 
earlier  panel  methods  by  employing  a  “higher  order”  panel  method;  that  is,  the  singularity  strengths  are  not 
constant  on  each  panel.  PANAIR  has  the  ability  to  analyse  arbitrary  geometry,  and  this  ability,  along  with 
the  ease  with  which  models  could  be  created  and  results  read  through  input  and  output  files,  greatly 
influenced  the  decision  to  select  it  as  the  flow  solver  to  utilise  in  the  simulations  and  coupling  with  the 
optimiser.  To  simplify  the  construction  of  PANAIR  input  files,  MAKEWGS  converts  the  set  of  data  points 
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describing  the  external  shape  of  the  wing  into  the  Langley  Wire  Frame  Geometry  Standard  (WGS) 
(Craidon  1985)  and  a  second  program,  PANIN,  combines  the  MAKEWGS  file  with  an  auxiliary  file  to 
generate  the  complete  PANAIR  input  file.  The  flow  conditions  contained  in  the  auxiliary  file  are  the 
altitude,  angle  of  attack  and  reference  data  for  the  wing. 

The  output  from  PANAIR  contains  the  calculated  pressures  and  corresponding  forces  and  moments  about 
the  simulated  geometry.  The  outputted  second  order  pressure  coefficients  were  selected  as  they  are  the 
preferred  ones  for  further  analysis  (Saaris  1992). 

PANAIR  seeks  to  solve  the  potential  flow  equation.  This  equation  is  a  combination  of  the  continuity 
equation,  momentum  and  speed  of  sound  equations,  subject  to  certain  conditions  such  as  the  flow  being 
irrotational.  After  considerable  mathematical  rearrangement,  these  equations  can  be  written  as: 
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This  equation  can  be  written  in  terms  of  the  perturbation  velocities,  (Eqn  2)  .The  right  hand  side  term  is 
often  neglected  as  it  is  of  second  order  in  the  perturbation  components. 
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By  bringing  in  the  velocity  potential  we  have: 
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PANAIR  solves  these  equations  and  uses  a  fully  continuous  quadratic  doublet,  and  linear  source 
distribution  along  both  the  span  and  chord-wise  directions. 

Complete  derivation  of  Eqn  4  can  be  found  in  Bertin  (Bertin  2002). 

In  PANAIR,  for  a  panel  solution  to  be  found  about  an  arbitrary  configuration  the  wing  or  similar  object 
must  be  closed,  i.e.  there  are  no  gaps  in  the  model.  This  means  that  for  the  simulations  performed,  the  wing 
described  has  to  be  mirrored  about  the  XZ  plane  such  that  it  made  a  complete  'flying  wing'.  Although  this 
model  description  would  not  yield  the  same  results  as  if  the  wing  were  attached  to  a  fuselage,  it  was 
concluded  that  for  ease  of  simulation,  this  method  would  be  chosen. 
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Figure  12  shows  the  upper  surface  of  the  aerodynamic  model  of  the  wing.  For  all  the  simulations 
performed,  the  authors  modelled  the  wings  reflected  about  the  XZ  plane  with  a  full  outer  skin  along  the 
upper  and  lower  surfaces  of  the  wing,  and  a  straight  wing  tip.  The  wing  tip  panels  are  plotted  in  Figure  13. 


A  wake  is  automatically  attached  to  the  trailing  edge  of  the  wing  by  PANAIR. 
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5.2.2  Validation  of  PANAIR 


Before  coupling  PANAIR  with  the  optimiser,  the  program  was  validated  against  experimental  wind  tunnel 
data  and  results  published  by  Margason,  et  al  (Margason,  Kjelgaard  et  al.).  The  test  case  considers  the  45° 
swept  back  wing  operating  at  M  ,  =  0.6  and  at  an  angle  of  attack  a  =  6.3  deg. 

Along  the  leading  and  trailing  edges,  the  panel  sizes  increased  exponentially  for  the  first  ten  percent  of  the 
chord,  and  then  remained  at  a  constant  size  for  the  remaining  eighty  percent.  Through  experimentation,  this 
was  found  to  yield  the  best  result  for  the  lowest  computational  cost. 

Results 

The  CPU  time  for  this  aerodynamic  computation  was  12  minutes  on  a  single  Pentium  4  2.4  GHz  processor. 


The  graded  coefficients  of  pressure  distribution  for  the  upper  and  lower  surfaces  are  shown  in  Figures  14 
and  15 
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Figure  15:  Pressure  Coefficients  about  the  Upper  Surface  of  the  Test  Section 


Figures  16  and  17  show  die  second  order  pressure  coefficient  distributions  at  two  span- wise  stations  as 
calculated  by  PANAIR.  The  results  from  the  fifty  five  percent  span  locations  are  shown  in  Figure  16.  As 
indicated  in  this  figure,  the  results  attained  are  in  close  agreement  with  the  experimental  results  reported  in 
References  (Kolbe  and  Boltz  1951)  and  (Margason,  Kjelgaard  et  ah). 


Figure  17  shows  the  results  at  8  %  span.  The  discrepancy  between  experimental  and  PANAIR  results  close 
to  the  wing  root  is  due  to  a  mismatch  in  the  Kutta  condition.  This  is  due  to  the  set  up  and  execution  of 
PANAIR  where  the  Kutta  condition  is  automatically  imposed  on  the  boundary  with  the  trailing  wing  wake. 
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This  resulted  in  an  incorrect  Kutta  condition  at  the  root  chord  and  hence  the  condition  of  the  upper 
pressures  equalling  the  lower  pressures,  as  should  be  found  on  the  leading  and  trailing  edges,  was  not 
correctly  realised.  Further  along  the  wing,  the  Kutta  condition  is  met  as  can  be  seen  from  Figure  1 6. 
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Concluding  the  validation  study,  it  is  shown  that  the  results  obtained  by  PANAIR  are  in  good  agreement 
with  experimental  data.  PANAIR  has  capabilities  to  provide  accurate  results  for  different  arbitrary 
geometries  and  solve  the  aerodynamic  characteristics  on  3D  wings.  PANAIR  provides  some  advantages 
over  more  complex  Navier-Stokes  solvers: 

•  good  accuracy  even  considering  the  linearilised  flow  assumption 

•  a  reduction  in  computational  expense. 

The  computational  time  for  a  single  PANAIR  calculation  varied  between  two  and  twelve  minutes  for  an 
arbitrary  wing  on  a  Pentium  4,  2.4  GHz  machine.  The  time  taken  was  dependant  on  the  number  of 
discritised  surface  panels  used  to  model  each  wing.  In  contrast,  a  full  Navier  Stokes  solution  takes  in  the 
order  of  several  hours.  Therefore,  the  aerodynamic  solver  has  the  capability  to  accurately  model  the  flow 
about  arbitrary  geometries  and  is  included  in  the  aero  structural  algorithm  and  overall  program  coupling 
with  the  optimiser. 

5.2.3  Friction 

As  mentioned  in  previous  section,  the  program  FRICTION  is  used  to  estimate  the  Form  and  Friction  drag 
produced  by  the  wing  at  the  simulation  conditions. 

Form  drag  is  drag  due  to  the  geometry  of  the  wing  being  analysed  and  is  influenced  by  the  frontal  area  of 
the  wing  ( t/c  ratio). 

Friction  drag  is  drag  due  to  the  movement  of  air  molecules  along  the  wing  surface  which  lead  to  the 
creation  of  boundary  layers.  If  the  boundary  layer  is  laminar,  momentum  and  energy  are  only  mixed 
between  neighbouring  layers  in  the  flow  and  on  a  microscopic  level;  hence  the  friction  drag  is  small.  As  the 
boundary  layer  becomes  turbulent  and  the  mixing  between  sections  of  the  boundary  layer  occurs  on  a 
macroscopic  level,  the  friction  drag  can  be  many  times  larger  than  that  for  a  laminar  boundary  layer. 
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FRICTION  calculates  the  laminar  boundary  layer  friction  coefficient  through  the  use  of  the  standard 
Blasius  formula  (Eqn  5)  with  the  inclusion  of  the  Chapman-Rubesin  constant  C  . 
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The  turbulent  boundary  layer  is  calculated  through  the  use  of  the  van  Driest  II  method.  For  a  full  details  on 
the  method  employed  when  calculating  the  turbulent  boundary  layer,  refer  to  Bertin  (Bertin  2002). 

Assumptions 

In  the  simulation  process  it  is  assumed  that  fifty  percent  of  the  wing  was  turbulent,  and  the  remainder 
laminar.  This  was  set  as  an  assumption  and  kept  for  all  simulations  as  smooth  graphite  epoxy  has  a  very 
low  roughness  value,  even  though  this  value  would  vary  between  simulated  wings  in  real  life.  The  value  of 
fifty  percent  chord  was  selected  after  consulting  pictures  of  ribbon  experiments  performed  on  numerous 
aircraft  wings  (Bertin  2002). 

5.3  Structural  Analysis 

The  structural  analysis  is  conducted  using  MSC.Nastran®,  a  modular  Finite  Element  Analysis  (FEA) 
software  package  written  by  MSC. Software"'.  MSC.Nastran"  has  modules  for  heat  transfer,  dynamics,  spot 
welding,  aero-elasticity,  and  nonlinear  analyses. 

In  this  work,  the  operation  of  MSC.Nastran  "  was  done  through  the  use  of  the  Bulk  Data  File  (BDF)  input 
method.  For  the  simulations  SOL  106  was  used  which  defines  the  Nonlinear  Static  analysis.  This  method 
was  chosen  over  SOL  101,  Linear  Static  as  it  was  expected  that  some  of  the  wing  geometries  defined  would 
undergo  large  deflections  and  or  large  strains.  It  was  therefore  necessary  to  correctly  determine  structural 
forces  formed  and  hence  the  Nonlinear  Static  solution  method,  as  described  in  Section  5.3.3,  was  chosen. 

5.3.1  Structural  Model  and  Constraints 


For  the  structural  analysis,  a  simplified  finite  element  model  constructed  from  a  varying  number  of  ribs  and 
spars  was  used.  The  majority  of  the  structural  model  was  assembled  using  shell  elements  with  the  spar  and 
rib  caps  modelled  as  rod  elements.  The  number  of  nodes  and  elements  and  hence  degrees  of  freedom 
(DOF)  varied  depending  on  the  internal  wing  structure.  A  sample  finite  element  model  is  illustrated  in 
Figure  18  and  19. 
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Figure  18:  Structural  finite  element  model  used  in  MSC.Nastran 
(Upper  surface  panels  -  blue.  Lower  surface  panels-  red) 

The  different  components  making  up  each  candidate  wing,  namely  spars,  ribs,  wing  skin  and  spar  and  rib 
caps,  are  described  in  more  detail  below. 

Spars 

The  Spars  were  linearly  distributed  over  the  chord  depending  on  their  number.  This  greatly  reduced  the 
number  of  variables  in  the  simulation,  but  did  not  necessarily  yield  the  optimal  result  for  any  given  number 
of  spars. 

The  thicknesses  of  the  Spars  were  modelled  to  decrease  according  to  a  parabolic  function  as  a  function  of 
the  span.  This  modelling  method  allowed  the  optimizer  the  ability  to  define  thicknesses  that  varied  at  a 
constant  rate  down  the  span  to  Spar  thickness  which  varied  parabolically  increasing  the  scope  of  the 
optimization  solution. 

Each  spar  was  broken  down  into  sections.  The  sections  were  defined  as  the  distance  between  consecutive 
ribs  and  this  allowed  for  a  minimum  of  MSC.Nastran8  staictural  elements  to  be  used  when  defining  the 
overall  structure  as  shown  in  Figure  19. 
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As  per  recommendations  by  Marisarla  (Marisarla  2005),  the  number  of  spar  element  sections  never 
decreased  below  8.  This  allowed  for  the  correct  modelling  of  nonlinear  element  bending  with  sufficient 
accuracy. 

Ribs 

As  with  the  Spars,  the  Ribs  were  also  broken  down  into  sections.  The  sections  were  defined  from  the 
leading  edge  to  the  first  spar,  between  consecutive  spars  if  they  existed  and  then  from  the  trailing  edge  spar 
to  the  wing  trailing  edge. 

The  Rib  thickness  was  also  varied  according  to  a  parabolic  function  in  a  similar  manner  as  to  the  Spar 
thicknesses. 

As  the  rounded  shape  of  the  leading  edge  provides  additional  torsional  support  in  the  structure,  the  leading 
edge  section  of  each  rib  was  modelled  using  a  quadrilateral  element  as  in  Figure  20.  Furthermore,  the 
trailing  edge  section  was  modelled  using  a  triangular  element  to  retain  the  correct  structural  modelling  as 
shown  in  Figure  21. 
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Figure  20:  Blunt  leading  Edge  Quadrilateral  Elements 


Wing  skin 

The  wing  skin  was  defined  by  quadrilateral  element  sections.  These  sections  covered  both  the  upper  and 
lower  surface  of  the  wing  and  were  primarily  used  as  a  means  of  accurately  placing  the  pressure  loads  onto 
the  internal  wing  structure.  Each  section  therefore  extended  between  neighbouring  spars  and  ribs,  and  in  the 
case  of  the  blunt  wing  leading  edge,  between  rib  points  as  in  Figure  22. 
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Figure  22:  Wing  Leading  Edge  Skin 


As  the  wing  skin  thickness  varies  from  wing  root  to  wing  tip  and  from  the  wing  leading  edge  to  the  wing 
trailing  edge,  so  the  wing  skin  thickness  in  the  MSC.Nastran”  structural  model  had  the  ability  to  have  a 
varying  wing  skin  thickness  as  well.  This  variation  was  a  simple  linear  function  for  both  the  variation  along 
the  span  and  chord. 

Spar  and  Rib  Caps 

Spar  and  rib  caps  significantly  add  to  the  structural  rigidity  of  the  wing  and  hence  were  modelled  as  lumped 
thicknesses  at  the  junction  between  spars  and  ribs  with  the  wing  skin  as  rod  elements.  As  with  the  rib  and 
spar  thicknesses,  the  rib  and  spar  caps  thicknesses  had  the  ability  to  taper  from  the  root  to  the  tip. 

5.3.2  Input  File  from  Structural  Model 

MSC.Nastran”  works  using  a  Bulk  Data  input  file  which  defines  the  type  of  solver  and  structural  model. 
Some  details  and  considerations  for  a  typical  BDF  file  in  the  simulation  are: 

Bulk  Data 

Before  the  definition  of  any  of  the  Element  Properties,  etc,  a  number  of  parameters  were  set  up  for  each 
run.  The  most  important  were  the  ones  affecting  nodal  constraint  and  solution  definition. 

As  the  simulation  being  performed  was  a  nonlinear  static  one,  the  requirement  to  constrain  nodes  which 
underwent  large  rotations  was  relaxed. 

The  Nonlinear  Static  solution  method  was  selected.  MSC.Nastran”1  simulations  were  run  with  fifty  load 
steps.  This  allowed  the  program  ample  chance  to  find  a  simulation  solution  within  a  minimum  time  frame. 
If  a  larger  number  were  chosen,  MSC.Nastran”  would  most  probably  find  most  simulation  solutions,  but 
the  cost  would  be  in  the  increase  in  computational  time  required  to  perform  all  the  structural  simulations  at 
this  small  load  step  size. 

The  other  feature  of  interest  in  the  nonlinear  parameter  setup  was  the  convergence  criteria  used.  As  the 
solutions  being  found  were  for  coarse  structural  models  defining  rough  design  analysis,  the  convergence 
criteria  were  relaxed.  This  was  backed  up  through  referring  to  the  ‘Quick  Reference  Guide’  for 
MSC.Nastran®  (MSC.Software  2005). 
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Element  Properties 


The  Spars,  Ribs  and  Wing  Skins  were  modelled  using  shell  elements.  Considerable  research  was  conducted 
into  different  modelling  techniques  which  have  been  used  when  defining  the  above  mentioned  structural 
elements  (Marisarla  2005),  (Venter  and  Sobieszczanski-Sobieski  2004).  It  was  decided  that  shell  elements 
would  be  used  to  define  the  structural  members  through  the  PSHELL  command.  For  these  structural 
elements  this  would  not  only  allow  bending  and  shear  within  each  of  the  elements,  but  would  also  account 
for  the  coupling  properties  such  as  Poisson’s  Ratio  which  would  not  be  accounted  for  if  only  shear 
elements  were  used. 

The  drawback  of  using  rod  elements  to  describe  the  lumped  rib  and  spar  caps,  described  in  later  sections,  is 
that  rod  elements  have  a  constant  cross  sectional  area.  Furthermore,  each  rod  element  with  a  different  cross 
sectional  area  needs  to  have  a  separate  entry  describing  the  properties  attributed  to  that  element.  This  was 
done  through  the  PROD  input  deck. 

Quadrilateral  Elements 

CQUAD4  input  decks  were  used  to  define  skin  panels,  wing  spars  and  wing  ribs  which  had  varying 
thickness  at  all  four  corners  of  the  quadrilateral  elements.  This  allowed  for  better  structural  modelling 
compared  to  the  use  of  a  constant  thickness  throughout  all  the  defined  element  sections. 

Furthermore,  CQUAD4  elements  define  plane  strain  plate  elements  and  hence  this  feature  was  beneficial  to 
the  structural  simulation. 

One  example  of  a  CQUAD  element  definition  is  as  follows: 

CQUAD4  31  1  498  816  812  494 

0.00058  0.00044  0.0004  0.00052 

This  method  of  CQUAD  description  allowed  for  the  definition  of  edge  thicknesses  at  node  points  along  the 
element. 

For  a  full  breakdown  of  the  CQUAD4  and  subsequent  input  decks,  please  refer  to  the  MSC.Nastran® 
reference  manuals  available  from  MSC. Software®. 

Triangular  Elements 

CTRI3  elements  were  used  to  model  the  trailing  edge  elements  defining  the  wing  ribs.  As  with  the 
CQUAD4  entries,  the  use  of  CTRI3  elements  allowed  for  the  specification  of  different  element  thicknesses 
at  each  of  the  node  points.  This  allowed  for  refined  structural  calculations. 

An  example  of  a  CTRI3  element  definition  is  as  follows: 

CTRIA3  252  1  1379  3340  1909 

0.00132  0.00132  0.00132 

The  same  definition  type  was  used  with  the  CQUAD  elements  as  was  used  with  the  CTRI3  elements  in  that 
it  allowed  the  authors  the  ability  to  define  the  edge  thicknesses  at  each  of  the  node  points. 

Rod  Elements 

CROD  elements  were  used  to  model  lumped  spar  and  rib  caps;  the  desirable  feature  being  that  they  are 
only  able  to  carry  tension  and  compression  loads  and  hence  do  not  alter  the  torsional  properties  of  the  wing 
section. 

An  example  of  a  CROD  element  definition  is  as  follows: 
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CROD, 266,4,266, 478 


This  entry  simply  defines  the  rod  element  according  to  the  nodes  located  at  each  end  of  the  element. 
Materials 

For  all  the  simulations  performed,  only  one  material  was  used.  Due  to  time  constraints  the  material  was 
defined  as  an  isotropic  graphite  epoxy  material  which  had  constant  material  properties  under  tension, 
compression  and  shear.  This  resulted  in  a  material  with  a  Poisson’s  Ratio  of  0.3,  Young’s  Modulus  of 
1.53E11  and  a  density  of  1310  kg/m3.  The  material  was  defined  through  the  MAT1  input  deck  and 
referenced  by  all  elements. 

An  example  of  MAT1  describing  the  Graphite  Epoxy  composite  used  with  both  the  High  and  Medium 
Altitude  Long  Endurance  vehicles  is  as  follows: 

MAT1 *  1  153000000000.0  0.30  *M1 

*M1  1310.01 

Due  to  the  large  values  used  when  describing  the  Young’s  Modulus  of  the  material,  the  long  field  format 
was  used. 

Nodes 

This  section  referenced  the  node  number  to  the  corresponding  x,  y  and  z  location.  The  node  numbers  were 
used  when  defining  the  CQUAD4,  CTRI3  and  CROD  elements. 

Loads 

Two  loads  were  applied  to  each  candidate  wing  structure.  The  second  order  pressure  coefficients  calculated 
by  PANAIR  and  converted  to  pressures  were  applied  to  the  upper  and  lower  wing  skin  panels,  and  a  gravity 
load  was  added  to  model  the  deflection  and  strains  due  to  the  mass  of  the  structure. 

The  gravity  load  was  required  as  the  mass  of  the  structure  would  create  a  downward  acting  force,  and  hence 
a  moment  about  the  X  axis,  which  would  then  counteract  the  lifting  force  and  moment  generated  by  the 
applied  pressures.  Therefore  the  gravity  load  would  decrease  the  deflection  of  the  wing  and  give  a  more 
realistic  indication  of  the  initial  deflection  under  the  applied  aerodynamic  loads. 

The  pressure  load  was  defined  through  the  use  of  the  PLOAD2  input  deck.  This  input  deck  defined  the 
magnitude  of  the  pressures  acting  perpendicular  to  the  surface  defining  a  wing  skin  panel.  Due  to  some 
pressures  being  larger  than  eight  characters  in  length;  the  authors  chose  to  use  the  large  field  format. 

PLOAD2*  2  -2917.363  54  *P54 

*P54 

From  the  above  example  code,  one  can  see  that  the  long  field  format  was  used  as  the  applied  pressures  on 
each  of  the  wing  skin  panels  in  column  three  could  become  very  large  about  the  leading  edges. 

Node  Constraints 

For  simplicity  it  was  decided  that  all  the  nodes  along  the  wing  root  would  be  constrained.  Although  this 
would  not  happen  in  the  case  of  a  wing  installed  into  the  fuselage  of  an  aircraft,  for  simulation  purposes 
and  to  ease  the  comparisons  made  between  results,  the  entire  wing  root  was  fixed  in  all  displacements  and 
rotations. 
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If  the  fuselage  was  included  in  the  simulation  environment,  the  wing  spars  would  simply  be  extended  to  the 
fuselage  centre  line  and  the  nodes  defining  the  ends  of  the  spars  fixed  in  both  displacements  and  rotations. 
This  would  be  the  correct  method  of  simulating  the  node  fixities  for  the  wing  as  the  fuselage  would  counter 
some  of  the  forces  and  moments  generated  by  the  wing. 

5.3.3  Use  of  MSC.Nastran®  for  Analysing  High  Aspect  Ratio  Wings 

As  it  was  unknown  to  what  amount  each  candidate  wing  would  deflect  under  the  calculated  aerodynamic 
loading,  it  was  necessary  to  use  a  generic  and  robust  method  when  solving  for  wing  deflection. 
Furthermore,  the  deflection  solution  would  need  to  be  accurate  if  it  were  to  be  included  in  determining  the 
overall  Pareto  front  members  and  if  the  variables  used  to  define  the  wing  were  to  be  passed  to  a  higher 
level  in  the  optimiser  structure.  It  was  therefore  decided  that  a  linear  static  solution  (Eqn  6),  calculated  in 
MSC.Nastran"  would  not  provide  accurate  results  as  this  method  would  only  solve  the  candidate  wing 
using  small  deflection  theory.  The  choice  of  solution  method  was  also  made  after  conclusions  drawn  by 
Marisarla  (Marisarla  2005). 


KQ=F  (6) 

As  some  of  the  candidate  wings  could  deflect  an  amount  greater  than  a  linear  theory  could  accurately 
model,  the  non-linear  solution  method  was  chosen.  The  limit  to  linear  theory  is  roughly  taken  in  the  case  of 
a  cantilever  beam  to  be  deflections  of  the  same  order  of  magnitude  as  the  beam  thickness.  The  change  to  a 
nonlinear  solution  method  only  required  changes  to  the  MSC.Nastran"  input  file. 

For  the  Nonlinear  Static  solution,  the  applied  loads  are  broken  down  into  a  number  of  sub  loads.  These 
loads  are  then  applied  in  steps  and  the  solution  energy  and  deflection  iterated  until  convergence  for  that 
load  step,  and  the  solution  sequence  continued  until  the  entire  load  has  been  applied  to  the  structure. 
MSC.Nastran"  has  the  ability  to  bisect  the  current  load  step  if  it  is  unable  to  find  a  solution  with  the  current 
step  size.  If,  even  after  a  number  of  user  specified  load  step  bisections,  MSC.Nastran"'  is  unable  to  find  a 
solution,  the  program  terminates  with  a  failure  message. 

Further  information  on  the  nonlinear  solution  method  used  by  MSC.Nastran  can  be  found  in  (Komzsik 

2001). 

A  draw  back  to  using  this  method  was  that  since  the  load  was  incrementally  increased  and  iterated  for 
deflection  and  energy  balance,  the  solution  time  was  longer  increasing  the  computational  cost  per  iteration. 

5.4  Aero-Structural  Analysis  of  Baseline  Designs 

To  illustrate  the  use  of  the  aero-  structural  analysis  we  apply  the  method  to  two  test  cases  related  to  UAV 
wing  design;  a  Medium  Altitude  Long  Endurance  (MALE)  research  UAV  similar  to  the  Altair  UAV  wing 
and  a  High  Altitude  Long  Endurance  (HALE)  UAVs  similar  to  the  Global  Hawk  wing. 

5.4.1  Medium  Altitude  Long  Endurance  (Male)  UAV 

Problem  Formulation 

A  detailed  analysis  of  a  wing  for  a  medium  altitude  research  UAV  application  similar  to  the  Altair  UAV 
was  considered.  The  operating  conditions  and  data  are  based  on  details  specified  by  NASA  (NASA  2004). 
The  aircraft  maximum  gross  weight  is  approximately  3175  kg,  has  a  wingspan  of  approximately  12.542  m, 
a  mean  chord  of  approximately  1.244  m,  and  a  planform  shape  with  2.25°  sweep  between  the  root  and 
second  crank  (Figure  23  and  Table  1). 
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Figure  23:  Design  Variables  for  Altair  (NASA  2004) 


An  assumption  was  made  that  the  aircraft  operated  at  210  KIAS  (108.10  m/s)  at  a  cruise  altitude  of 
41,600//. 

It  was  assumed  that  the  wing  uses  the  NACA4415  aerofoil  at  all  spar  locations  as  information  on  the  actual 
aerofoils  used  was  unavailable.  Table  2  summarises  some  of  the  flight  properties  for  the  aircraft.  These 
conditions  assume  an  aircraft  at  mid  weight-cmise  during  an  extended  cruise  phase  at  intermediate  altitude. 


Characteristics,  Units 

Maximum  Altitude,  (m) 

15850 

Cruise  Altitude  (80%MaxAlti),  (nt) 

12680 

Density  (©,  Max  Altitude,  (kg/m3) 

0.16492 

Density  (h  Cruise  Altitude,  (kg/m3) 

0.27830 

Air  Speed,  (m/s) 

108.10 

Speed  of  Sound  @  Max  Altitude,  (m/s) 

295.073 

Speed  of  Sound  (a/Cruisc,  (m/s) 

295.073 

Mach  @Max  Altitude 

Mm  =  0.3663 

Mach  @Cruise 

Moo  =  0.3663 

Gross  Weight,  (kg) 

3175 

Table  1:  Flight  Conditions  for  Altair 


Minimum  Required  Lift  Coefficient 

The  following  equation  was  used  to  calculate  the  required  coefficient  of  lift  at  a  number  of  flight  conditions 
(Eqn  7). 
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Wx  9.81 


(7) 


C‘=T 


,2  x  p  x  V~  x  S 

The  results  after  using  Eqn  7  are  summarised  in  Table  2. 


Flight  Condition 

Required  CL 

Max  Gross  Take  Off  at  Maximum  Altitude 

1.1066 

Beginning  Cruise  after  T/O  and  Climb  (5%  fuel  consumption) 

0.6417 

Middle  of  Cruise  (50%  fuel  left) 

0.5152 

right  before  Landing  (10%  fuel  left  ->  reserve) 

0.4028 

Table  2:  Required  Altair  CL  at  different  Flight  Conditions 


Simulations  were  performed  at  an  angle  of  attack  of  4°  which  yielded  a  calculated  Cl  greater  than  0.64. 
This  is  equivalent  to  operating  the  Altair  UAV  at  the  beginning  of  the  cmise  condition  as  tabulated  in  Table 
2. 

Aerodynamics  Model 

The  PANAIR  and  FRICTION  codes  were  used  for  the  aerodynamic  analysis  and  are  described  in  Section 
5.2.  In  total,  1224  aerodynamic  panels  were  used  to  discretise  the  Altair  MALE  model.  Figure  24  shows  the 
external  planform  shape. 
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Figure  24:  External  Planform  layout  of  Altair  UAV 


Figure  25  shows  the  linear  (mid  chord  section)  and  exponential  (leading  and  trailing  edges)  layout  of  the 
quadrilateral  panels  comprising  the  upper  surface  of  the  aerodynamics  model. 
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Structural  Model 

The  structural  model  is  described  by  the  values  displayed  in  Table  3  and  shown  in  Figure  26. 


Variable 

Value 

Number  of  spars 

2 

Number  of  ribs 

10 

Materials 

Wing  skin 

Graphite/epoxy 

Spars 

Graphite/epoxy 

Ribs 

Graphite/epoxy 

Thicknesses 

Wing  skin  (m) 

0.0001 

Spars (m) 

0.08 

Ribs  (m) 

0.005 

Thickness  Ratios 

Wing  Skin  (Tip) 

0.1 

Wing  Skin  (Trailing  Edge) 

0.1 

Spars 

0.05 

Ribs 

0.25 

Areas 

Spar  Cap  (m2) 

0.00375 

Rib  Cap  (m2) 

0.0015 

Material  Properties 

Youngs  Modulus  (Pa) 

1.53E11 

Poissons  Ratio 

0.3 

Density  (A#/m3) 

1310 

Table  3:  Structural  Variable  Values  for  Altair  UAV 
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Results 


Aerodynamic  analysis 

Figures  27  and  28  show  the  wing  span  pressure  coefficient  distribution.  The  overall  Altair  wing 
aerodynamic  characteristics  for  the  above  defined  configuration  are  listed  in  Table  4. 


The  wing  was  found  to  weigh  220.9  kg  which  is  in  close  agreement  with  the  value  estimated  from  Raymer 
(Raymer  1999). 


L/D 

23.9 

CL 

0.66478 

cD 

0.028 

cm 

-0.1768 

Deflection 

12.06% 

Table  4:  Altair  Wing  Aerodynamic  and  Deflection 
Characteristics 


Figure  27  shows  the  pressure  distribution  over  the  upper  wing  surface  and  Figure  28  that  of  the  lower 
surface. 
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Figure  27:  Pressure  distribution  along  Altair  Wing  Top  surface 
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Figure  28:  Pressure  distribution  along  Altair  Wing  Bottom  surface, 
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Performing  a  Mach  sweep  where  the  angle  of  attack  for  the  vehicle  was  kept  constant  yielded  the  results  in 
Figure  29  for  the  calculated  induced  component  of  the  total  CD  value. 
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Figure  29:  CD  results  for  Altair  Mach  Sweep  between  0.2  and  0.4 


Performing  the  same  calculations  for  the  coefficient  of  lift  yielded  the  results  shown  in  Figure  30. 
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Figure  30:  Cl  results  for  Altair  Mach  Sweep  between  0.2  and  0.4 


Structural  Analysis 

Figures  31  and  32  show  the  displacement  and  strain  results  respectively  after  applying  the  pressure  values 
calculated  using  PANAIR. 
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MSC.Patran  2005  1 0-Dec+05 1 7:51 :40 
Fringe:  DEFAULT.SC1.  A1  :Non-linear:  1 00.  %  of  Load,  Displacements,  Translational,  Magniti 
Deform:  DEFAULT.SC1,  A1:Non-linear:  100.  %  of  Load,  Displacements,  Translational, 


1 .57+000 
1  5 1 .47+000 

ION-^Y§g|^o 
1 .26+000 
1.15+000| 

1 .05+0001 
9.44-001  T 
8.39-001  [ 
7.34-001 1 
6.29-001 1 
5.24-001 
4.19-001 
3.15-OOU 

2.1 0-001 E 

1.05-001 1 

0.1 

default_Fringe : 

Max  1 .57+000  @Nd  665 
Min  0.  @Nd96 
default_Deformation : 
Max  1 .57+000  @Nd  665 

Figure  31:  Non-linear  Displacement  for  Altair  UAV  wing  under  2.5g  load  case 


MSC.Patran  2005  1 0-Deo05  1 7:59:34 

Fringe:  DEFAULT.SC1,  ATNon-linear:  100.  %  of  Load.  Nonlinear  Strains,  Strain  Tensor, .  At  Center 
Deform:  DEFAULT.SC1.  A1:Non-linear:  100.  %  of  Load,  Displacements.  Translational, 


4.03-004 

3.76-004 


default_Fringe  : 

Max  4.03-004  @Nd  898 
Min  0.  @Nd  96 
default_Deformation : 


Max  1 .57+000  @Nd  665 


Figure  32:  Non-linear  Strain  for  Altair  UAV  wing  under  2.5g  load  case 


The  calculations  above  resulted  in  the  definition  of  the  benchmark  for  the  Altair  MALE  UAV  redesign  to 
be  considered  in  Section  7. 1 . 
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5.4.2  High  Altitude  Long  Endurance  (Hale)  UAV  Wing 


Problem  Formulation 

As  with  the  Altair  MAVE  UAV  a  detailed  analysis  was  performed  on  a  UAV  similar  to  the  Global  Hawk 
HALE.  The  operating  conditions  and  data  are  based  on  the  information  provided  in  Air  Force  Technology 
(www.airforce-technology.com  2004).  The  aircraft  maximum  gross  weight  is  approximately  12,100  kg,  has 
a  wingspan  of  approximately  17.7  m,  a  mean  chord  of  approximately  1.72  m,  and  a  planform  shape  with 
5.9°  sweep  between  the  root  and  tip  (Figure  33).  The  assumption  was  made  that  the  aircraft  operated  at  334 
KIAS  (176.56  m/s)  at  a  cruise  altitude  of  15850  m. 

From  Selig  (Selig  2002),  the  Global  Hawk  makes  use  of  the  NASA  LRN  1015  aerofoil  section  throughout 
the  wing. 


Figure  33:  Planform  layout  of  Global  Hawk  (Goraj,  Frydrychewicz  et  al.  2004) 


The  Global  Hawk  flight  conditions  are  summarised  in  Table  5. 


Characteristics,  Units 

Maximum  Altitude,  (m) 

19812 

Cruise  Altitude  (80%MaxAlti),  (m) 

15850 

Density  (a).  Max  Altitude,  (kg/m3) 

0.08761 

Density  (©Cruise  Altitude,  (kg/m3) 

0.16938 

Air  Speed,  (m/s) 

176.566 

Speed  of  Sound  (a  Max  Altitude,  (m/s) 

295.073 

Speed  of  Sound  @Cruise,  (m/s) 

295.073 

Mach  @Max  Altitude 

Mo„  =  0.5983 

Mach  (©Cruise 

Moo  =  0.5983 

Gross  Weight,  (kg) 

11,600 

Table  5:  Flight  Conditions  for  the  Global  Hawk 
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Minimum  Required  Lift  Coefficient 

Eqn  7  was  used  to  calculate  the  required  coefficient  of  lift  at  each  flight  condition  indicated  in  Table  5.  The 
results  after  using  Eqn  7  are  summarised  in  Table  6. 


Flight  Condition 

Required  CL 

Max  Gross  Take  Off  at  Maximum  Altitude 

1.7341 

Beginning  Cruise  after  T/O  and  Climb  (5%  fuel  consumption) 

0.8962 

Middle  of  Cmise  (50%  fuel  left) 

0.6718 

right  before  Landing  (10%  fuel  left  ->  reserve) 

0.4723 

Table  6:  Required  Altair  Cl  at  different  Flight  Conditions 


Simulations  were  performed  at  an  angle  of  attack  of  4.75°  which  yielded  a  calculated  Cl  greater  than 
0.8962.  This  is  equivalent  to  operating  the  Global  Hawk  UAV  at  the  beginning  of  the  cmise  condition  as 
tabulated  in  Table  6. 

Aerodynamics  Model 

The  PANAIR  and  FRICTION  codes  were  used  for  the  aerodynamic  analysis  and  are  described  in  Section 
5.2.  In  total,  1980  aerodynamic  panels  were  used  to  discritise  the  Global  Hawk  HALE  benchmark  model. 
Figure  34  shows  the  linear  (mid  chord  section)  and  exponential  (leading  and  trailing  edges)  layout  of  the 
quadrilateral  panels  comprising  a  portion  of  the  upper  surface  of  the  aerodynamics  model.  This  is  the  same 
model  layout  as  was  used  with  the  Altair  benchmark  in  Section  5.4.1. 


Structural  Model 

The  structural  model  is  described  by  the  values  displayed  in  Table  7  and  shown  in  Figure  36.  Figure  35 
shows  the  internal  structure  of  the  wing  and  Figure  36  the  external  panels. 
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Variable 

Value 

Number  of  spars 

5 

Number  of  ribs 

10 

Materials 

Wing  skin 

Graphite/epoxy 

Spars 

Graphite/epoxy 

Ribs 

Graphite/epoxy 

Thicknesses 

Wing  skin  (m) 

0.01 

Spars  (m) 

0.08 

Ribs  (m) 

0.0015 

Thickness  Ratios 

Wing  Skin  (Tip) 

0.01 

Wing  Skin  (Trailing  Edge) 

0.01 

Spars 

0.01 

Ribs 

0.01 

Areas 

Spar  Cap  (m2) 

0.01 

Rib  Cap  (m2) 

0.0005 

Material  Properties 

Youngs  Modulus  (Pa) 

1.53E11 

Poissons  Ratio 

0.3 

Density  (kg/m3) 

1310 

Table  7:  Structural  Variable  Values  for  Global  Hawk  UAV 
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Figure  36:  External  Planform  layout  of  Global  Hawk  UAV 


Results 

Aerodynamic  analysis 

The  overall  Global  Hawk  wing  aerodynamic  characteristics  for  the  above  defined  configuration  are  listed  in 
Table  8.  Figures  37  and  38  show  the  wing  span  pressure  coefficient  distribution. 


L/D 

31.729 

CL 

0.905 

cD 

0.02849 

cm 

-0.208 

Deflection 

16.7% 

Table  8:  Global  Hawk  Wing  Aerodynamic  and 
Deflection  Characteristics 


The  wing  was  found  to  weigh  1138.15  kg  which  is  in  close  agreement  with  the  value  calculated  through  the 
use  of  Raymer  (Raymer  1999). 

Figure  37  shows  the  pressure  distribution  over  the  upper  wing  surface  and  Figure  38  that  of  the  lower 
surface. 
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Figure  37:  Pressure  distribution  along  Global  Hawk  Wing  Top  surface 
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Figure  38:  Pressure  distribution  along  Global  Hawk  Wing  Bottom  surface 
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Performing  a  Mach  sweep  where  the  angle  of  attack  for  the  vehicle  was  kept  constant  yielded  the  results  in 
Figure  39  for  the  calculated  induced  component  of  the  total  CD  value. 


0.022 


0.021 


Q  0.02 


0.019 


■&  0.018 
o 


0.017 


0.016  I - 

0.45 


0.5 


i_ i 

0.55  0.6 

Mach  Number 


0.65 


0.7 


Figure  39:  CD  results  for  Global  Hawk  Mach  Sweep  between  0.475  and 


Performing  the  same  calculations  for  the  coefficient  of  lift  yielded  the  results  shown  in  Figure  40. 


Mach  Number 

Figure  40:  CL  results  for  Global  Hawk  Mach  Sweep  between  0.475  and 


Structural  Analysis 

Figures  41  and  42  show  the  displacement  and  strain  results  respectively  after  applying  the  pressure  values 
calculated  using  PANAIR. 
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Figure  41:  Non-linear  Displacement  for  Global  Hawk  UAV  wing  under  2.5g  load  case. 
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Figure  42:  Non-linear  Strain  for  Global  Hawk  UAV  wing  under  2.5g  load  case 


The  calculations  above  resulted  in  the  definition  of  the  benchmark  for  the  Global  Hawk  HALE  UAV 
redesign  to  be  considered  in  Section  7.2 


51 


6  Aero-Structural  Optimisation 


In  a  general  form,  the  aero-structural  optimisation  process  follows  the  flow  diagram  in  Figure  43.  The 
Aero-Structural  Solver  computes  the  flow  dynamics  for  the  planform  and  structural  variables  defined  by 
the  Evolutionary  Optimiser  and  create  the  four  output  files  described  in  Section  5.1.2. 


Figure  43:  Aero-structural  optimisation  flow  chart 


In  Figure  43,  the  Evolutionary  Optimiser  uses  a  hierarchical  approach  with  three  levels,  on  the  bottom  level 
a  coarse  analysis  is  performed  to  direct  the  exploration,  at  the  top  level  a  more  precise  model,  that  better 
describes  the  physics  involved,  is  used  and  at  an  intermediate  level,  a  compromised  balance  between  the 
top  and  bottom  layers  is  used.  Initially  the  system  will  specify  the  design  variables,  constraints  and 
parameters  and  will  generate  a  random  sub  population  of  individuals  at  each  layer.  It  then  defines  the 
number  of  sub  populations  (nodes)  and  number  of  hierarchical  levels  which  for  simplicity  is  equal  to  the 
number  of  analyses  being  performed.  Once  these  initial  populations  are  generated,  the  algorithm  will 
initialise  all  the  required  populations  through  repeated  calls  to  the  Aero- Structural  Solver  and  go  through  a 
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series  of  steps.  The  scheduler  first  determines  whether  the  given  stopping  conditions  have  been  met,  and  if 
so,  the  evolutionary  loop  is  exited  and  the  entire  process  is  stopped.  If  no  stopping  conditions  are  met,  the 
scheduler  updates  the  asynchronous  solver  so  that  further  progress  can  be  made.  The  scheduler  determines 
whether  or  not  candidate  solutions  which  have  been  solved  are  ready  for  incorporation  into  the  population. 
If  such  solutions  exist,  the  incorporation  routine  is  called  and  available  candidates,  which  now  have  had  a 
fitness  assigned  (evaluated  by  the  Aero- Structural  Solver),  are  processed;  it  receives  the  individual,  ages 
the  population  and  buffer,  performs  Pareto  tournament  selection,  deletes  the  oldest  from  the  buffer  and  if 
the  acceptance  is  true,  it  is  inserted  in  the  population  which  was  subsequently  sorted.  It  then  updates  the 
CMA  parameters.  Finally  the  scheduler  determines  whether  it  is  possible  to  generate  more  candidates  by 
polling  the  asynchronous  solver.  If  this  is  possible,  then  the  generation  routine  is  called  and  individuals  are 
generated  via  the  evolutionary  operators  through  recombination  and  mutation  via  CMA  always  checking 
the  upper  and  lower  bounds.  During  evaluation,  the  optimiser  will  take  output  analyses  and  parameters  to 
guarantee  the  satisfaction  of  constraints  and  compute  the  overall  fitness  function.  If  the  problem  is  multi 
objective,  the  algorithm  will  find  the  non-dominated  individuals  and  will  calculate  the  Pareto  fronts. 

As  MSC.Nastran8  runs  on  the  availability  of  licenses  on  a  master  server,  the  Aero-Structural  Solver  was 
written  such  that  the  code  checks  to  see  whether  or  not  a  license  file  exists  when  attempting  to  analyse  a 
candidate  wing  with  MSC.Nastran8  .  If  a  license  does  not  exist  as  all  the  licenses  are  currently  being  used 
by  other  analyses,  the  program  waits  for  a  license  to  become  available  and  does  not  simply  discard  the 
simulation. 

This  ability  to  wait  for  licenses,  therefore  allows  the  Evolutionary  Optimiser  to  be  run  on  more  machines 
than  there  are  MSC.Nastran1'  licenses.  This  is  beneficial  as  only  a  fraction  of  the  total  time  taken  to  solve 
each  candidate  wing  solution  is  taken  calculating  the  structural  solution.  Therefore,  if  only  the  same 
numbers  of  computers  as  MSC.Nastran8  licenses  are  run,  there  will  be  instances  when  MSC.Nastran8 
licenses  are  idle.  This  ability  to  wait  for  licenses  therefore  is  beneficial. 

6.1  Handling  of  Constraints 

Constraints  are  handled  using  penalty  functions  which  are  broken  up  into  two  different  parts.  The  first  part 
penalized  the  aerodynamics  of  the  wing,  and  the  second  the  structure.  It  was  found  to  be  beneficial  to 
separate  the  penalties  as  it  meant  that  if  the  optimizer  generates  a  candidate  solution  which  resulted  in  a 
light  and  rigid  structure,  while  the  aerodynamics  were  far  from  optimum,  the  structure  could  still  exist 
within  the  Pareto  set  and  have  its  variable  values  passed  onto  further  generations  of  candidate  solutions. 
The  same  consideration  applied  to  the  aerodynamics  of  a  candidate  wing.  If  the  aerodynamics  were  found 
to  be  near  optimal  while  the  wing  deflection  large  with  failed  structural  components  due  to  high  strain 
values,  it  would  slow  down  the  optimizer  if  the  solution  were  merely  discarded. 

Therefore,  by  programming  the  penalty  function  to  only  penalize  those  components  which  failed,  the 
overall  efficiency  of  the  optimizer  was  increased. 

To  further  assist  the  EA  when  evaluating  potential  non-dominated  solutions,  an  exponential  penalty 
function  is  used.  This  forces  penalised  wings  into  regions  outside  the  Pareto  front,  or  off  to  one  extreme 
favouring  only  one  of  the  fitness  functions.  This  meant  that  only  wings  which  were  not  penalised,  or 
penalised  slightly,  existed  close  to  the  Pareto  front.  This  would  not  be  seen  with  the  use  of  a  linear  penalty 
function  as  a  wing  which  was  penalised  for  defecting  50%  could  be  seen  as  a  wing  with  the  same  weight  as 
a  200%  wing. 

6.1.1  Aerodynamics 

Aerodynamic  moment 

In  order  to  keep  the  control  system  of  the  vehicle  relatively  unchanged,  the  authors  decided  to  impose 
penalties  on  the  tested  wing  if  the  aerodynamic  moment  generated  by  the  wing  was  greater  than  the  original 
benchmark  tests  performed. 
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The  wing  was  linearly  penalized  by  increasing  the  fitness  value  defined  by  the  inverse  of  the  lift  to  drag 
ratio  for  every  percent  over  the  aerodynamic  moment  benchmark  value  calculated. 

Coefficient  of  Lift 

If  the  optimizer  were  to  design  a  wing  with  aerodynamics  such  that  the  inverse  of  the  lift  to  drag  ratio  was 
minimized,  and  only  the  aerodynamic  moment  penalized,  the  code  might  evolve  a  wing  which  would 
produce  too  little  lift  when  placed  in  the  required  flight  regime.  Furthermore,  if  one  were  then  to  alter  the 
angle  of  attack  for  the  aerofoil  to  attain  the  required  Coefficient  of  Lift,  the  total  drag  for  the  aerofoil  would 
change  and  hence  there  lies  a  possibility  that  the  aerofoil  would  cease  to  be  the  optimal  configuration. 

The  authors  therefore  decided  to  penalize  a  candidate  wing  if  the  lift  created  was  below  a  set  amount.  This 
amount  was  usually  set  at  ninety  percent  of  the  required  lift  at  the  start  of  the  cruise  phase  of  the  aircraft. 
As  with  the  Aerodynamic  moment  penalty,  the  Coefficient  of  Lift  penalty  was  added  to  the  fitness  function 
defined  by  the  inverse  of  the  lift  to  drag  ratio  at  a  rate  for  every  percent  under  the  required  value.  If  the 
Coefficient  of  Lift  calculated  for  the  wing  was  above  the  required  value,  no  penalty  was  imposed. 

6.1.2  Structure 

Aerofoil  Thickness 

The  thickness  of  an  aerofoil  shape  has  implications  not  just  for  the  aerodynamics  of  the  vehicle.  If  a  very 
thin  aerofoil  section  is  designed  it  might  produce  excellent  aerodynamic  results  when  built  and  tested,  but 
there  would  be  a  very  large  drawback  in  that  the  space  between  the  upper  and  lower  surfaces  of  the  wing 
would  be  very  small  resulting  in  a  cramped  and  tight  working  environment.  This  small  working 
environment  would  affect  the  internal  layout  of  the  wing  and  the  ease  of  component  inspection. 

Therefore  the  thickness  of  each  aerofoil  must  exceed  a  value  of  12%  ( t/c  >  0.12).  If  the  Evolutionary 
Optimiser  produced  a  candidate  aerofoil  which  had  a  smaller  thickness  value,  the  aerofoil  would  be 
penalised  by  equally  penalising  both  the  aerodynamic  and  structural  values  via  a  linear  penalty  method.  In 
addition,  aerofoil  sections  generated  by  the  Evolutionary  Optimiser  which  had  thicknesses  outside  some 
bounds,  usually  between  10%  and  15%,  were  rejected  immediately  before  any  analysis  was  performed  with 
them. 

Buckling  of  Wing  Skin  Panels 

Thin  sheets  of  material  under  compression  and  or  shear  loads  are  subject  to  buckling.  When  a  sheet  of 
material  buckles  under  the  applied  load,  the  magnitude  of  the  strength  with  which  the  sheet  of  material  is 
able  to  withstand  reduces  significantly. 

An  accurate  method  of  predicting  whether  or  not  any  of  the  elements  within  a  candidate  wing  would  be 
subject  to  buckling  would  be  to  perform  a  buckling  analysis  in  MSC.Nastram  .  The  additional 
computational  overhead  would  be  rather  large  increasing  the  time  required  to  determine  the  fitness  of  each 
candidate  wing.  This  would  further  extend  the  time  required  by  HAPMOEA  to  find  a  global  optimal 
solution  given  the  simulation  bounds.  It  therefore  became  appropriate  to  find  another  method  to  determine 
whether  or  not  any  elements  within  a  candidate  wing  would  buckle.  As  the  Spars  and  Ribs  were  of 
considerable  thicknesses  compared  to  the  wing  skin,  and  also  since  the  loads  and  moments  carried  by  these 
two  structural  elements  lay  in  plane  with  the  elements  as  in  Figure  44a,  the  possibility  of  any  of  the 
sections  making  up  the  Spars  and  Ribs  buckling  would  be  very  small  compared  to  the  possibility  of  the  skin 
panel  sections  buckling.  This  is  due  to  the  fact  that  the  applied  loads  and  moments  across  the  Skin  sections 
lie  out  of  plane  as  in  Figure  44 b. 
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Following  the  method  followed  by  Venter  and  Sobieszczanski-Sobieski  (Venter  and  Sobieszczanski- 
Sobieski  2004),  simple  analytical  expressions  were  used  to  test  whether  or  not  the  stresses  in  the  Wing  Skin 
sections  were  large  enough  to  cause  buckling.  These  local  buckling  expressions  are  described  by  Eqn.  8. 


Where: 


1<0 

°cr 

^L-1<0 


(  \ 


-  + 


-1<0 


cr  cr  J 


N, 


3  ME 


2 1 ,C7cr  fb 


v 1  y 

Nn  4.85  kE 


It 


b 

vO 


(8) 


and  b  is  the  width  of  each  Wing  Skin  section  and  k  was  an  additional  safety  factor  to  guard  against 
delamination  of  the  plies  making  up  the  carbon  graphite  composite  material.  The  safety  factor  was  set  at 
four  for  all  the  optimisations. 

Tip  Twist 

As  the  wing  bends,  the  pressures  generated  by  the  wing  would  change  and  hence  the  forces  within  the  wing 
would  also  change.  If  the  wing  tip  twisted,  there  is  a  possibility  that  the  wing  would  greatly  increase  the 
drag  created  by  the  wing,  and  also  increase  the  moment  about  the  Z  axis.  As  a  full  aero-elastic  analysis  was 
not  part  of  the  simulation  runs,  any  further  increase  in  the  angle  at  which  the  wing  would  make  with  the 
freestream  flow  would  be  disadvantageous  to  the  computation  as  the  final  wing  position  and  internal  forces 
could  not  be  computed. 
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It  was  therefore  necessary  to  penalize  wings  in  which  after  the  structural  simulation  using  MSC.Nastran8, 
the  wing  tip  angle  of  attack  was  greater  than  a  user  set  value  from  the  initial  wing  setting  angle.  This 
penalty  took  the  form  of  mass  added  to  the  structure  at  the  rate  for  each  additional  degree. 

Deflection 

As  no  full  aeroelastic  analysis  was  performed  on  each  candidate  wing,  the  final  aerodynamic  and  structural 
solution  was  not  known.  Therefore  the  assumption  was  made  that  if  the  initial  deflection  of  the  wing  was 
small,  the  wing  would  have  a  smaller  chance  of  aerodynamic  divergence  if  the  aerodynamic  and  structural 
simulations  were  performed  on  the  deflected  wing  was  made.  Although  this  might  not  be  the  case  for 
certain  wings,  it  was  simple  and  easy  to  implement  and  placed  a  far  smaller  overhead  on  the  computational 
cost  compared  to  the  analysis  of  the  mass  matrix  for  each  candidate  wing. 

The  wing  was  therefore  allowed  to  deflect  a  set  percentage  of  the  span  for  a  2.5 g  load  case  before  a  penalty 
value  was  added  to  the  fitness  value.  As  with  the  Tip  Twist  condition,  the  penalty  took  the  form  of 
additional  mass  and  was  added  per  percent  deflection  over  twenty  percent. 

Failure  of  Internal  Component  Panels 

In  this  aero-structural  analysis  the  testing  of  structural  members  for  failure  was  performed  through  the  use 
of  the  maximum  strain  theory.  This  was  done  as  stress  varies  between  plies  and  orientations  in  a  composite 
material,  while  strain  varies  linearly  through  the  thickness.  This  same  method  was  employed  in  Venter  and 
Sobieszczanski-Sobieski  (Venter  and  Sobieszczanski-Sobieski  2004). 

If  any  components  making  up  the  internal  structure  within  the  wing,  spars,  ribs,  skin,  etc  failed  due  to 
excessive  strains,  the  wing  was  heavily  penalized  by  the  addition  of  extra  wing  mass.  This  mass  was 
increased  at  a  rate  of  ten  percent  for  each  panel  making  up  a  component  which  failed.  Hence,  if  ten 
CQUAD4  components  making  up  a  spar  failed,  the  mass  of  the  wing  structure  was  increased  one  hundred 
percent. 

6.1.3  Combined  Penalties 

If  a  wing  was  analysed  and  the  penalties  were  very  large  against  a  user  set  maximum,  instead  of  only 
penalising  one  of  the  fitness  functions,  both  fitness  functions  would  be  penalised  according  to  the 
maximum  penalty  associated  with  either  of  the  calculated  values.  This  forces  the  solution  to  be  dominated 
and  hence  the  genetic  information  is  not  passed  on  to  any  further  generations  of  candidate  wing  solutions. 
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7  Aero-Structural  Optimisation  Test-Cases 


This  section  considers  the  application  of  the  optimisation  methodology  for  two  test  cases  related  to  UAV 
design.  The  first  test  case  is  for  a  Medium  Altitude  Long  Endurance  (MALE)  UAV  and  the  second  for  a 
High  Altitude  Long  Endurance  (HALE)  UAV. 

7. 1  Medium  Altitude  Long  Endurance  (MALE)  UAV  Design  And 

Optimisation 

This  first  test  case  considers  the  multidisciplinary,  multi-objective  design  optimisation  of  a  MALE  UAV 
Wing  similar  to  the  Altair  MALE  Unmanned  Aerial  Vehicle  (UAV)  considered  in  Section  5.4.1.  The  two 
objectives  are  the  maximisation  of  the  lift  to  drag  ratio  (L/D),  and  the  minimisation  of  wing  weight.  The 
cruise  Mach  number  and  altitude  are  0.3663  and  12680  m  respectively  and  the  wing  area  is  set  to  29.21  nr. 

7.1.1  Design  Variables  and  Constraints. 

The  wing  geometry  is  represented  by  three  aerofoil  sections  with  sixteen  variables  for  each  section,  eight 
variables  that  describe  the  wing  planform  and  eleven  variables  to  describe  the  internal  structure. 

The  aerofoil  geometry  is  represented  by  the  combination  of  a  mean  line  and  thickness  distribution,  which  is 
a  very  common  procedure  in  classical  aerodynamics  (Abott  and  Doenhoff  1980).  Both  lines  are  represented 
by  Bezier  curves  with  leading  and  trailing  edge  points  fixed  at  (0.0, 0.0)  and  (1.0, 0.0)  respectively,  and  a 
variable  number  of  intermediate  control  points  whose  x-positions  are  fixed  in  advance  and  whose  v-heights 
form  the  problem  unknowns.  In  this  case  6  free  control  points  on  the  mean  line  and  10  free  control  points 
on  the  thickness  distribution  are  taken. 

The  wing  planform  design  variables  are  indicated  in  Figures  45  and  46  and  their  upper  and  lower  bounds 
are  represented  in  Table  9.  In  total  67  design  variables  are  used  for  the  MALE  UAV  wing  optimisation. 
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Description 

Benchmark 

Lower  Bound 

Upper  Bound 

Angle  of  Attack,  a  (deg) 

4 

0 

6 

Wing  Semi  Span,  b  (m) 

12.542 

8.2 

18.5 

Root  Chord,  cr  (m) 

1.646 

1.6 

1.7 

Wing  Inboard  Leading  Edge  Sweep,  Arc  (rad) 

0.039 

0 

5 

Wing  Outboard  Leading  Edge  Sweep,  Act  (rad) 

0.039 

0 

5 

Chord  Ratio  at  Break,  Xrc 

0.444 

0.35 

0.65 

Chord  Ratio  at  Tip,  Xct 

0.444 

0.35 

0.65 

Crank  Location,  crank 7 

0.731 

0.5 

0.8 

Number  of  Spars,  Ns 

2 

1 

3 

Number  of  Ribs.  Nr 

10 

8 

12 

Rib  Root  Thickness,  Rrt  (m) 

0.005 

0.003 

0.0065 

Rib  Root  Thickness  Taper  Ratio,  RrTr 

0.25 

0.15 

0.3 

Spar  Root  Thickness,  Srt  (m) 

0.08 

0.05 

0.09 

Spar  Root  Thickness  Taper  Ratio,  SrTr 

0.05 

0.01 

0.2 

Wing  Skin  Root  Thickness,  WsRt  (m) 

0.001 

0.0005 

0.002 

Wing  Skin  Thickness  Tip  Taper  Ratio,  WstTr 

0.1 

0.05 

0.15 

Wing  Skin  Thickness  Trailing  Taper  Ratio,  WsTre 

0.1 

0.05 

0.15 

Spar  Cap  Root  Area,  Sc  (nr) 

0.00375 

0.0025 

0.04 

Rib  Cap  Root  Area,  Rc(m2) 

0.0015 

0.001 

0.002 

Table  9:  Structural  and  Planform  Design  Variables  for  the  MALE  UAV  wing 
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7.1.2  Design  Constraints  and  Penalties 


There  are  several  constraints  which  must  be  satisfied  for  each  candidate  wing  namely: 

•  The  pitching  moment  must  not  be  greater  than  -0.1768. 

•  The  calculated  coefficient  of  lift  must  be  greater  than  0.64. 

•  The  thickness  of  each  aerofoil  must  exceed  12%  (t/c  >  0. 12). 

•  No  internal  or  external  structural  sections  can  buckle  or  fail  due  to  excessive  strains. 

•  The  wing  tip  may  not  twist  more  than  one  degree. 

•  The  wing  may  not  deflect  more  than  twenty  percent  of  the  span. 

All  the  above  mentioned  constraints  are  applied  by  penalising  either,  or  both  fitness  values  via  an 
exponential  penalty  method  as  described  in  Section  6.1.  In  addition,  any  aerofoil  generated  outside  the 
thickness  bounds  of  10%  to  15%  are  rejected  immediately  before  any  analysis  is  performed  on  the  wing. 
The  above  data  is  summarised  in  Table  10. 


Description 

Values 

Allowable  Strain 

0.00333 

Allowed  wing  tip  twist  in  degrees 

1 

Allowed  wing  deflection  as  a  percentage  span 

20 

Allowed  wing  moment  (benchmark) 

-0.1768 

Minimum  lift  to  be  generated  by  the  wing 

0.64 

Wing  Mass  per  degree  twist  of  wing  tip 

0.1 

Wing  Mass  per  deflection  percent  over  20 

0.1 

Wing  Mass  per  failed  panel 

0.1 

Additional  CD  per  over  allowable  aerodynamic  moment 

0.001 

Additional  CD  per  less  than  the  required  minimum  Cl 

0.005 

Penalise  both  fitness  functions  above  a  Drag  ratio  of 

0.1 

Penalise  both  fitness  functions  above  a  Mass  ratio  of 

0.1 

Table  10:  Constraint  and  Penalty  Values 


7.1.3  Fitness  Functions 

The  two  fitness  functions  to  be  optimised  are  defined  as  Eqn  9  and  Eqn  10: 

min  (/, ) :  fx=  +  Penalty  ^ 

min  (/2) :  f2=W  +  penalty  (10) 

Eqn  9  relates  to  the  aerodynamics  of  the  candidate  wing.  The  Evolutionary  Optimiser  seeks  to  minimise  the 
inverse  of  the  lift  to  drag  ratio  with  the  addition  of  penalties  which  are  described  in  Section  6. 1 . 

Eqn  10  relates  to  the  structural  aspects  of  the  candidate  wing.  This  function  has  the  units  of  mass  and  seeks 
to  minimise  the  calculated  structural  mass  of  the  wing.  Penalties  are  also  added  to  Eqn  10  which  are 
described  in  Section  6.1 

For  both  Eqns  9  and  10,  the  penalties  are  calculated  based  on  the  constraints  specified  in  Section  7.1.2. 

7.1.4  Aero-structural  Analysis 

The  aerodynamic  and  structural  characteristics  of  the  wing  configurations  are  evaluated  using  the  Aero- 
Structural  Solver  described  in  Section  5.1. 
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7.1.5  Implementation 


The  solution  to  this  problem  has  been  calculated  using  two  approaches;  the  first  one  uses  a  traditional  EA 
method  with  a  single  population  model  based  upon  a  fine  meshed  aerodynamic  model  of  each  candidate 
wing.  Between  1904  and  2992  aerodynamic  panels  were  used  to  define  the  candidate  wings  using  a  fine 
mesh  structure.  The  second  method  uses  a  hierarchical  topology  of  mesh  resolutions  comprising  three 
levels.  The  same  aerodynamic  resolution  as  was  used  in  the  single  population  experiments  was  used  for  the 
top  level,  a  more  coarse  mesh  model  used  in  an  intermediate  aerodynamic  level  (between  1428  and  2244 
aerodynamic  panels),  and  a  coarse  mesh  (between  952  and  1496  aerodynamic  panels)  used  in  the  lowest 
model  to  explore  the  design  search  space. 

This  resulted  in  the  coarse  mesh  solutions  requiring  roughly  4  minutes  to  solve  for  the  aerodynamic 
pressures  about  a  candidate  wing,  and  15  minutes  for  a  fine  mesh.  And  intermediate  wing  took  roughly  10 
minutes  on  a  single  2.4  GHz  machine. 

Three  P4,  2.4  GHz  machines  were  used  in  the  calculation  and  the  population  size  for  this  problem  was  set 
to  30  in  all  levels.  This  was  set  for  both  the  single  and  hierarchical  solutions. 

For  both  solution  approaches,  the  single  population  and  the  hierarchical  method,  the  structural  model  panel 
density  was  kept  as  a  function  of  the  internal  geometry  and  was  not  based  upon  the  level  of  model 
complexity.  For  a  description  of  the  method  used,  please  refer  to  Section  5.3. 

7.1.6  Results 

The  algorithm  was  run  for  500  function  evaluations  and  took  36  hours  to  compute  on  a  cluster  of  three 
machines.  Figure  47  shows  the  Pareto  fronts  obtained  by  using  the  two  approaches.  The  Pareto  members 
selected  for  further  investigation  are  also  indicated.  By  inspecting  Figure  47  it  can  be  seen  that  the  use  of  a 
multi-fidelity  (hierarchical)  approach  gives  an  overall  lower  front  as  compared  to  a  single  model  approach. 
Both  approaches  give  results  which  lie  to  the  left  of  the  benchmark  value  for  the  MALE  wing. 
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Figure  47:  Pareto  Fronts  for  MALE  UAV  Optimisation 
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The  algorithm  progress  for  the  hierarchical  approach  is  shown  in  Figure  48.  The  fitness  for  all  levels 
decreases  with  time  and  it  can  be  seen  that  after  350  evaluations,  the  fitness  of  each  level  reaches  a  steady 
state  value  indicating  an  optimum  solution  given  the  constrains  imposed. 


0.027 


0.026 


0.025 


0.024 


0.023 


0.022 


0.021 


0.02 


! - { - ! - ! - I - F~ 


"-4-.  ’•! 


!  ! _ 1! 

i 

\ — 

“T  — 


■  I 
1 


T' 

i  1 


—  Level  1 

-  Level  2 

-  Level  2 

—  Level  3 
■■  Level  3 
■■  Level  3 
■■  Level  3 


50  100  150  200  250  300  350  400  450  500  550 

Function  Evaluations 

Figure  48:  MALE  Hierarchical  Algorithm  Progress 


61 


The  optimisation  was  run  for  a  second  time  but  with  only  one  population  to  test  the  validity  of  the 
hierarchical  procedure  employed  by  the  HAPMOEA  algorithm.  The  results  are  shown  in  Figure  49. 

The  algorithm  progress  shown  in  Figure  49  relates  to  the  Pareto  front  shown  as  a  solid  red  line  on  Figure  47. 


Figure  49  MALE  Single  Population  Algorithm  Progress 


From  Figures  47  to  49  it  can  be  seen  that  the  performance  for  the  Hierarchical  scheme  was  higher  than  that 
of  the  single  population  scheme  over  the  same  number  of  function  evaluations  of  the  Level  1  node.  The 
hierarchical  approach  therefore  yields  better  solutions  for  the  same  number  of  function  evaluations  of  the 
Level  1  node. 

Three  wings  are  taken  for  further  evaluation  (Pareto  members  3,  8  and  23)  from  the  hierarchical 
optimisation  strategy  results  to  illustrate  the  two  objective  extremes  and  a  compromise  between  the  two 
results.  According  to  the  optimiser,  Pareto  3  has  the  best  aerodynamics,  Pareto  23  the  best  structural 
response  to  the  aerodynamic  loading  and  8  lies  between  the  two  solutions. 

Table  1 1  summarises  and  compares  the  values  for  the  design  variables  and  the  objective  functions  for  these 
three  geometries  against  the  benchmark. 
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Description 

Benchmark 

PM  3 

PM  8 

PM  23 

Angle  of  Attack,  a  (deg) 

4 

4.055 

3.957 

5.16798 

Wing  Semi  Span,  b  (m) 

12.542 

13.871 

13.385 

12.4694 

Root  Chord,  cr  (m) 

1.646 

1.638 

1.641 

1.66618 

Wing  Inboard  Leading  Edge  Sweep,  Arc  (rad) 

0.039 

0.0871 

0.0694 

0.0801 

Wing  Outboard  Leading  Edge  Sweep,  Ac,  (rad) 

0.039 

0.0255 

0.0486 

0.0105 

Chord  Ratio  at  Break,  Xrc 

0.444 

0.422 

0.448 

0.427 

Chord  Ratio  at  Tip, 

0.444 

0.475 

0.388 

0.477 

Crank  Location,  crank / 

0.731 

0.583 

0.667 

0.787 

Number  of  Spars,  Ns 

2 

1 

2 

2 

Number  of  Ribs.  Ns 

10 

11 

10 

9 

Rib  Root  Thickness,  Rrt  (m) 

0.005 

0.00351 

0.00465 

0.00435 

Rib  Root  Thickness  Taper  Ratio,  RrTr 

0.25 

0.182 

0.287 

0.243 

Spar  Root  Thickness,  Srt  (m) 

0.08 

0.0805 

0.0576 

0.0571 

Spar  Root  Thickness  Taper  Ratio,  SrTr 

0.05 

0.131 

0.154 

0.019 

Wing  Skin  Root  Thickness,  WsRt  (m) 

0.001 

0.000885 

0.000953 

0.00199 

Wing  Skin  Thickness  Tip  Taper  Ratio,  WstTr 

0.1 

0.122 

0.1 

0.0541 

Wing  Skin  Thickness  Trailing  Taper  Ratio,  WsTre 

0.1 

0.106 

0.0608 

0.108 

Spar  Cap  Root  Area,  Sc  (m2) 

0.00375 

0.00908 

0.00314 

0.00314 

Rib  Cap  Root  Area,  Rc(m2) 

0.0015 

0.00123 

0.00175 

0.00105 

Table  11:  Summary  and  comparison  of  design  variables  for  MALE  UAV 


Figure  50  shows  a  top  view  of  the  wings  selected  for  further  analysis  and  Figures  51  to  53  compare  the 
baseline  aerofoils  and  the  aerofoil  sections  of  each  Pareto  member  at  the  root,  crank  and  tip. 
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Figure  50:  Top  Views  of  Pareto  members  selected  for  further  analysis 
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On  inspecting  Figure  50,  the  optimiser  found  that  the  original  layout  of  the  wing  favoured  the  structural 
constraints  more  than  the  aerodynamics  and  changed  the  wing  span  and  crank  location  between  the 
different  Pareto  members  to  achieve  better  aerodynamics. 


X 

Figure  51:  Selected  MALE  Pareto  member  Root  Aerofoil  Sections 


X 

Figure  52:  Selected  MALE  Pareto  member  Break  Aerofoil  Sections 
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It  can  be  seen  that  classical  aerodynamic  shapes  have  been  evolved,  even  considering  that  the  optimisation 
was  started  completely  from  random  and  the  evolution  algorithm  had  no  problem  specific  knowledge  of 
appropriate  solution  types. 

Of  interest  is  that  all  the  aerofoil  sections  generated  by  the  HAPMOEA  code  and  shown  in  Figures  5 1  to  53 
have  a  lower  t/c  ratio  than  that  of  the  benchmark  aerofoil  section.  There  are  only  slight  differences  between 
the  root  aerofoil  sections  describing  the  Pareto  members  selected  in  Figure  51. 

A  Mach  number  sweep  was  performed  at  a  constant  angle  of  attack  for  the  selected  Pareto  members  and  the 
results  are  displayed  in  Figure  54  for  coefficient  of  drag  and  Figure  55  for  coefficient  of  lift.  The 
benchmark  wing  is  included  for  comparison. 
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Figure  54:  Coefficient  of  Drag  vs.  Mach  number  for  selected  MALE 
Pareto  members  and  Benchmark 


Over  the  range  of  Mach  numbers  investigated,  the  Pareto  members  have  lower  drag  coefficient  than  the 
benchmark  wing.  The  Pareto  members  three  and  eight  operate  at  an  almost  constant  drag  value  over  the 
range  of  Mach  numbers. 


As  the  Lift  penalty  function  described  in  Section  6.1.1  only  penalised  candidate  wing  which  produced  a  lift 
coefficient  lower  than  the  required  value  to  maintain  the  cruise  altitude  and  velocity,  the  selected  Pareto 
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member  wings  all  produce  more  than  the  required  lift.  This  extra  lift  value  is  small  through  compared  to  the 
overall  required  force,  being  roughly  0.02  for  Pareto  member  eight  compared  to  the  benchmark. 


Coefficient  of  Lift 

Figure  56:  Coefficient  of  Lift  vs.  Coefficient  of  Drag  for  selected 
MALE  Pareto  members  and  Benchmark. 


When  the  Coefficient  of  Lift  is  plotted  versus  the  Coefficient  of  Drag  for  the  selected  Pareto  members  and 
the  benchmark  MALE  wing  in  Figure  56,  one  can  see  that  the  performance  of  all  the  wings  belonging  to 
the  Pareto  front  is  higher  than  that  of  the  benchmark  wing.  The  required  Coefficient  of  Lift  of  0.64  is  met  at 
a  lower  Coefficient  of  Drag  with  all  the  selected  Pareto  member  wings. 

The  selected  MALE  UAV  Pareto  member  are  further  analysed  by  performing  an  Angle  of  Attack  sweep 
between  three  and  six  degrees  and  displayed  in  Figures  57  and  58.  For  all  the  analyses  performed,  the  Mach 
number  was  kept  constant  at  0.3663,  the  operating  velocity  of  the  benchmark  vehicle. 
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For  all  the  wings  tested,  the  Coefficient  of  Lift  shown  in  Figure  57  varied  in  almost  a  linear  fashion  over 
the  angles  tested.  Furthermore,  all  the  wings  behaved  in  a  similar  fashion  in  that  the  gradients  of  the  Lift 
slopes  are  all  very  similar. 
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As  with  the  Mach  sweep  performed  and  displayed  in  Figure  54,  the  selected  Pareto  members  displayed  and 
tested  in  Figure  58  have  a  lower  Coefficient  of  Drag  over  all  tested  angles  of  attack.  Between  the  angles  of 
three  and  four,  Pareto  23  and  the  bench  exhibit  the  same  reduction  in  drag.  Pareto  eight  does  not  exhibit  a 
similar  drop  over  the  angle  range  covered  while  Pareto  three  exhibits  the  reduction  between  four  and  four 
and  a  half  degrees. 

Figures  59  and  60  plot  the  second  order  pressure  coefficients  for  the  selected  Pareto  wings. 
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Figure  59:  External  Pressure  Coefficients  for  selected  MALE  Pareto 
Members  -  Upper  Surface 


69 


Pareto  23 


Figure  60:  External  Pressure  Coefficients  for  selected  MALE  Pareto 
Members  -  Lower  Surface 


In  comparison  with  the  MALE  UAV  benchmark  wing  pressure  coefficients  displayed  in  Figures  27  and  28, 
the  wings  selected  from  the  Pareto  set  produce  higher  pressures  along  the  leading  edge  of  the  wings.  The 
maximum  coefficient  of  pressure  for  the  benchmark  wing  is  negative  three  while  the  Pareto  members  all 
produce  coefficients  of  pressure  as  high  as  negative  five.  This  could  be  due  to  the  smaller  t/c  ratio  of  the 
Pareto  aerofoils  evident  in  Figures  51  to  53.  The  pressures  generated  along  the  lower  surface  of  the  Pareto 
wings  are  also  smaller  in  magnitude  than  that  of  the  benchmark. 

A  larger  resultant  lift  force  is  therefore  created  by  the  Pareto  member  wings  and  this  is  in  good  agreement 
with  the  high  performance  of  the  wings  noted  when  the  Coefficients  of  Lift  and  Drag  were  plotted  in  Figure 
56. 

Figures  61  to  63  show  the  deflected  Pareto  wings.  The  original  wing  position  and  internal  structure  is 
shown  as  a  wire  frame  for  comparison. 
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Figure  61:  MALE  Pareto  3  Displacement 
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Figure  62:  MALE  Pareto  8  Displacement 
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Figure  63:  MALE  Pareto  23  Displacement 


All  the  MALE  wings  deflect  nonlinearly.  From  the  Figures  above  it  can  be  seen  that  the  deflection  at  the 
wing  tips  is  greater  than  the  thickness  of  the  aerofoil  sections  and  hence  wing  spars.  The  nonlinear  analysis 
solutions,  selected  for  the  MSC.Nastran8  structural  calculations,  was  therefore  necessary  as  the  deflections 
encountered  are  too  large  for  a  Linear  Static  analysis  and  if  used,  would  have  yielded  incorrect  results. 

An  analysis  of  the  strains  produced  by  the  aerodynamic  forces  is  shown  in  Figures  64  to  66. 
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Figure  64:  MALE  Pareto  3  Strain 
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Figure  65:  MALE  Pareto  8  Strain 


The  aerodynamic  load  generated  by  the  aerofoil  sections  cause  the  candidate  wings  to  want  to  rotate  about 
the  wing  root  (X  axis).  As  the  wings  were  fixed  in  all  directions  and  all  rotations  at  the  wing  root,  the 
maximum  strain  does  not  occur  at  the  wing  root  but  rather  a  distance  from  the  root,  where  the  rotation 
occurs. 
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Pareto  3 

Pareto  8 

Pareto  23 

Benchmark 

Deflection  (%) 

16 

17.3 

13.8 

12.5 

Maximum  Strain 

0.000396 

0.00052 

0.000485 

0.000403 

Mass  (kg) 

233.735 

217.873 

185.1 

220.9 

L/D 

43.367 

44.4 

44.25 

23.9 

Mass  Penalty  (kg) 

40 

0 

0 

0 

Drag  Penalty 

0.0027 

0 

0 

0 

Table  12:  Summary  of  Displacement,  Strain  and  Mass  for  Selected  MALE  Pareto  members 

and  the  Benchmark 


A  summary  of  each  of  the  Pareto  and  benchmark  wing  overall  characteristics  is  shown  in  Table  12  for 
comparison. 

It  can  be  noted  from  Table  12  that  all  the  Pareto  member  wings  out-performed  the  benchmark  wing  on  L/D 
ratio  with  Pareto  members  8  and  23  had  weights  which  were  lower  than  that  of  the  benchmark  wing  when 
the  penalty  values  were  included  in  the  overall  wing  mass.  All  the  wings  had  strain  values  which  were 
below  the  maximum  allowable  strain  defined  in  Table  10. 

7.2  High  Altitude  Long  Endurance  (Hale)  UAV  Wing  Design  And 
Optimisation 

In  the  second  case,  the  detailed  analysis  and  optimisation  of  a  wing  for  a  high  altitude  long  endurance 
(HALE)  UAV  application  is  considered.  This  is  a  multi-objective  problem  where  the  goal  is  to  maximise 
the  lift-to-dag  ratio  and  minimise  wing  weight.  The  operating  conditions  and  data  for  the  baseline  wing  are 
based  on  reference  (NORTHROP-GRUMMAN  2006)The  aircraft  has  a  wingspan  of  approximately  17  m,  a 
mean  chord  of  approximately  2.43  m,  a  wing  area  of  50.19  m "  and  a  plan  form  shape  with  5.9  deg  sweep. 
A  breakdown  of  the  benchmark  wing  can  be  found  in  Table  13. 

The  aircraft  is  assumed  to  operate  at  Mro  =  0.5983  at  a  cruise  altitude  of  15840  m.  It  is  also  assumed  that 
this  UAV  uses  a  single  aerofoil,  the  LRN1015  aerofoil,  as  the  only  wing  section  aerofoil  along  the  wing 
span. 

7.2.1  Design  Variables  and  Constraints. 

Similar  to  the  previous  test  case,  the  wing  geometry  is  represented  by  three  aerofoil  sections  consisting  of 
sixteen  variables  each,  eight  variables  that  describe  the  wing  planfonn  and  eleven  variables  to  describe  the 
internal  structure.  Once  again,  6  free  control  points  on  the  mean  line  and  10  free  control  points  on  the 
thickness  distribution  were  used  to  describe  each  of  the  aerofoil  sections. 

The  benchmark,  upper  and  lower  bounds  for  the  design  variables  are  indicated  in  Table  13. 
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Description 

Benchmark 

Lower  Bound 

Upper  Bound 

Angle  of  Attack,  a  (deg) 

5.75 

0 

6 

Root  Chord,  cr  (m) 

2.432 

2.375 

2.5 

Wing  Inboard  Leading  Edge  Sweep,  Arc  (rad) 

0.102972 

0 

0.1745 

Wing  Outboard  Leading  Edge  Sweep,  Act  (rad) 

0.102972 

0 

0.1745 

Chord  Ratio  at  Break, 

0.78 

0.55 

0.9 

Chord  Ratio  at  Tip,  lc, 

0.287 

0.2 

0.8 

Crank  Location,  crank / 

0.0881 

0.06 

0.6 

Number  of  Spars,  Ns 

5 

3 

7 

Number  of  Ribs.  Ns 

16 

12 

18 

Rib  Root  Thickness,  Rrt  (m) 

0.0015 

0.001 

0.005 

Rib  Root  Thickness  Taper  Ratio,  RrTr 

0.05 

0.01 

0.15 

Spar  Root  Thickness,  Srt  (m) 

0.09 

0.045 

0.091 

Spar  Root  Thickness  Taper  Ratio,  SrTr 

0.05 

0.01 

0.15 

Wing  Skin  Root  Thickness,  WsRt  (m) 

0.001 

0.0001 

0.05 

Wing  Skin  Thickness  Tip  Taper  Ratio,  WstTr 

0.01 

0.0095 

0.15 

Wing  Skin  Thickness  Trailing  Taper  Ratio,  WsTre 

0.01 

0.0095 

0.15 

Spar  Cap  Root  Area,  Sc  (m2) 

0.012 

0.004 

0.0125 

Rib  Cap  Root  Area,  Rc(m2) 

0.0005 

0.00045 

0.002 

Table  13:  Structural  and  Planform  Design  Variables  for  the  HALE  UAV  wing 


7.2.2  Design  Constraints  and  Penalties 

The  constraints  which  must  be  satisfied  for  each  candidate  wing  are  indicated  in  Table  14  along  with  the 
rates  with  which  wings  which  fail  the  constraints  are  penalised. 


Description 

Values  Allowed 

Allowable  Strain 

0.00333 

Allowed  wing  tip  twist  in  degrees 

1 

Allowed  wing  deflection  as  a  percentage  span 

20 

Allowed  wing  moment  (benchmark) 

-0.3041 

Minimum  lift  to  be  generated  by  the  wing 

0.89 

Wing  Mass  per  degree  twist  Penalty  Values 

0.1 

Wing  Mass  per  degree  over  20  span 

0.1 

Wing  Mass  per  failed  panel 

0.1 

Additional  Cd  per  over  allowable 

0.001 

Additional  Cd  per  less  than  the  required  minimum 

0.005 

Penalise  both  fitness  functions  above  a  Drag  ratio  of 

0.1 

Penalise  both  fitness  functions  above  a  Mass  ratio  of 

0.1 

Table  14:  Constraints  and  Penalty  Rates 


All  the  above  mentioned  constraints  are  applied  by  penalising  either  fitness  values  via  a  linear  penalty 
method  as  described  in  Section  6.1.  In  addition,  any  aerofoil  generated  outside  the  thickness  bounds  of 
10%  to  15%  are  rejected  immediately  before  any  analysis  is  performed  on  the  wing. 

7.2.3  Fitness  Functions 

The  two  fitness  functions  to  be  optimised  are  defined  as  Eqn  1 1  and  Eqn  12: 
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min  (/i): 

min(/2): 


f  =  (l  /  d)  + penalty 


f2=W  +  penalty 


(ID 

(12) 


Eqn  1 1  relates  to  the  aerodynamics  of  the  candidate  wing.  The  Evolutionary  Optimiser  seeks  to  minimise 
the  inverse  of  the  lift  to  drag  ratio  with  the  addition  of  penalties  which  are  described  in  Section  6.1. 

Eqn  12  relates  to  the  structural  aspects  of  the  candidate  wing.  This  function  has  the  units  of  mass  and  seeks 
to  minimise  the  calculated  mass  of  the  wing.  Penalties  are  also  added  to  Eqn  12  which  are  described  in 
Section  6.1 

For  both  Eqns  11  and  12,  the  penalties  are  calculated  based  on  the  constraints  specified  in  Section  7.1.2. 

7.2.4  Aero-structural  Analysis 

The  aerodynamic  and  structural  characteristics  of  the  wing  configurations  are  evaluated  using  the  Aero- 
Structural  Solver  described  in  Section  5.1  and  represented  in  Figure  67. 


7.2.5  Implementation 

A  single  optimisation  was  performed  for  the  HALE  UAV  wing.  This  optimisation  used  a  hierarchical 
approach  with  three  levels  and  seven  populations.  The  top  level  (Level  1)  used  a  fine  mesh  of  between 
2992  and  4624  aerodynamic  panels  to  calculate  the  pressures  about  the  candidate  wings.  A  more  coarse 
mesh  model  was  used  in  the  intermediate  aerodynamic  level  (between  2244  and  3468  aerodynamic  panels), 
and  a  coarse  mesh  (between  1496  and  2312  aerodynamic  panels)  used  in  the  lowest  mode  to  explore  the 
design  search  space. 
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Three  P4,  2.4  GHz  machines  were  used  in  the  calculation  and  the  population  size  for  this  problem  was  set 
to  30  in  all  levels  and  for  both  the  single  and  hierarchical  solutions. 

The  structural  model  panel  density  was  kept  as  a  function  of  the  internal  geometry  and  was  not  based  upon 
the  level  of  model  complexity. 

7.2.6  Results 

The  algorithm  was  run  for  300  function  evaluations  and  took  in  average  three  days  to  compute  on  a  cluster 
of  three,  2.4  GHz  machines.  Figure  68  compares  some  of  the  members  of  the  Pareto  front  and  the  baseline 
geometry.  Figure  69  shows  the  algorithm  progress  for  objective  one,  each  step  in  the  figure  roughly 
corresponds  to  migration  step-  good  information  from  the  bottom  levels  have  been  seeded  to  the  upper 
levels.  Figures  71  to  73  compare  the  aerofoils  at  root  break  and  tip  for  each  of  the  members  of  the  Pareto 
front  and  the  aerofoil  for  the  baseline  geometry.  Table  15  indicates  the  final  values  of  design  variables  and 
objective  functions  for  the  baseline  geometry  and  some  members  of  the  Pareto  front. 
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Figure  68:  Pareto  Selected  HALE  Members  and  Benchmark 


The  Pareto  selected  HALE  members  are  shown  against  the  benchmark  mass  and  inverse  of  the  L/D  ratio  in 
Figure  68.  From  the  above  Figure,  the  selected  Pareto  members  outperform  the  benchmark  wing. 

The  algorithm  progress  for  the  HALE  UAV  wing  optimisation  is  shown  in  Figure  69.  It  can  be  seen  that  the 
Level  1  fitness  decreases  with  function  evaluations  performed. 
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Function  Evaluations 

Figure  69:  Evolution  progress  and  migration  steps 


Table  15  summarises  and  compares  the  values  for  the  design  variables  for  the  three  selected  geometries 
against  the  benchmark. 


Description 

Benchmark 

PM  0 

PM  7 

PM  14 

Angle  of  Attack,  a  (deg) 

5.75 

5.427 

5.716 

5.428 

Wing  Semi  Span,  b  (m) 

16.941 

17.768 

16.586 

16.199 

Root  Chord,  cr  (m) 

2.432 

2.403 

2.401 

2.389 

Wing  Inboard  Leading  Edge  Sweep,  Arc  (rad) 

0.103 

0.111 

0.0769 

0.0996 

Wing  Outboard  Leading  Edge  Sweep,  Act  (rad) 

0.103 

0.0141 

0.0799 

0.0165 

Chord  Ratio  at  Break,  Xrc 

0.78 

0.744 

0.804 

0.845 

Chord  Ratio  at  Tip,  Xct 

0.287 

0.255 

0.283 

0.255 

Crank  Location,  crank i 

0.0881 

0.124 

0.117 

0.14 

Number  of  Spars,  Ns 

5 

4 

3 

3 

Number  of  Ribs.  Ns 

16 

14 

15 

15 

Rib  Root  Thickness,  Rrt  (m) 

0.0015 

0.00122 

0.0024 

0.0013 

Rib  Root  Thickness  Taper  Ratio,  RrTr 

0.05 

0.023 

0.0215 

0.03 

Spar  Root  Thickness,  Srt  (m) 

0.09 

0.0906 

0.091 

0.0904 

Spar  Root  Thickness  Taper  Ratio,  SrTr 

0.05 

0.128 

0.133 

0.12 

Wing  Skin  Root  Thickness,  WsRt  (m) 

0.001 

0.000975 

0.0103 

0.00119 

Wing  Skin  Thickness  Tip  Taper  Ratio,  WstTr 

0.01 

0.0645 

0.0962 

0.0942 

Wing  Skin  Thickness  Trailing  Taper  Ratio,  WsTre 

0.01 

0.0569 

0.0391 

0.0501 

Spar  Cap  Root  Area,  Sc  (m2) 

0.012 

0.008 

0.00655 

0.00508 

Rib  Cap  Root  Area,  Rc(m~) 

0.0005 

0.000867 

0.00131 

0.00102 

Table  15:  Summary  and  Comparison  of  HALE  UAV 
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Planform  Shapes 

Figure  70  shows  the  planform  shapes  for  the  selected  HALE  Pareto  members  along  with  the  Benchmark. 
From  the  Figure  one  is  able  to  see  that  all  the  wings  have  a  high  aspect  ratio  and  that  the  reference  area  was 
kept  constant. 

The  evolutionary  optimiser  produced  wings  which  seemed  to  favour  wings  with  a  lower  leading  edge 
sweep  angle.  The  crank  location  between  the  Pareto  members  and  the  Benchmark  wing  are  all  at  a  rather 
consistent  point  in  the  wing. 


Figure  70:  HALE  Pareto  and  Benchmark  planform  shapes 


Aerofoil  Sections 

The  evolutionary  optimiser  determined  a  number  of  different  aerofoil  sections  when  constructing  Pareto 
members.  The  root,  crank  and  tip  aerofoil  sections  are  shown  below  for  Pareto  members  0,  7  and  14  along 
with  the  Benchmark  LRN  1015  aerofoil  section. 
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Figure  72:  Comparison  of  Crank  aerofoil  geometries  for  selected  Pareto 
Members  and  baseline 
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From  the  above  Figures  71  to  73,  one  can  see  the  large  variation  in  aerofoil  shapes.  The  variations  though 
seem  to  have  a  common  trend  as  at  the  trailing  edges,  all  the  aerofoils  rise  above  the  leading  and  trailing 
edges  and  then  descent  to  the  trailing  edge.  The  high  variation  in  aerofoil  shapes  is  consistent  with  the 
variation  in  planform  shapes  shown  in  Figure  70  as  each  variation  requires  different  aerodynamic 
characteristics  in  an  attempt  to  meet  the  requirements  of  the  penalty  functions. 

Also  of  interest  as  with  the  sweep  angles  in  Figure  70,  the  Pareto  aerofoil  shapes  have  a  lower  upper 
surface  height  compared  to  the  benchmark  aerofoil. 

It  must  be  noted  that  the  Pareto  member  aerofoil  sections  all  have  a  lower  t/c  ratio  compared  to  that  of  the 
benchmark  aerofoil  shape. 

The  diverse  range  of  aerofoil  shapes  further  confirms  that  the  use  of  a  Pareto  front  when  evaluating  the 
HALE  UAV  wing  problem  was  a  success  at  there  are  large  variations  in  geometry. 

Mach  Sweep 

A  Mach  number  sweep  was  performed  for  the  selected  Pareto  members  and  the  results  are  displayed  in 
Figure  74  for  coefficient  of  drag  and  Figure  75  for  coefficient  of  lift. 
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Mach  Number 

Figure  74:  Comparison  of  CD  Mach  sweep  of  selected  HALE  Pareto 
members  and  benchmark 


The  Pareto  member  wings  produce  less  drag  than  the  benchmark  wing  over  the  range  of  Mach  numbers 
investigated.  Pareto  member  zero  produced  the  lowest  drag  at  the  operational  Mach  number  of  0.59,  of 


0.02175. 


Mach  Number 

Figure  75:  Comparison  of  CL  Mach  Sweep  of  selected  HALE  Pareto 
members  and  baseline  geometry 
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For  both  the  Coefficients  of  Lift  and  Drag,  the  trends  in  the  Figures  for  the  Pareto  members  are  very  similar. 


From  Figure  76,  the  performance  of  the  Pareto  selected  wings  can  be  compared  to  that  of  the  benchmark 
HALE  UAV  wing.  As  was  indicated  in  the  Pareto  front  shown  in  Figure  68,  Pareto  zero  has  the  best  L/D 
ratio  and  Pareto  fourteen  the  worst.  Even  though  Pareto  fourteen  has  the  best  structural  wing  solution,  it 
still  has  a  better  aerodynamic  performance  than  the  benchmark  wing. 

The  second  order  external  aerodynamic  pressures  are  plotted  for  the  upper  and  lower  Pareto  wing  surfaces 
in  Figures  77  and  78  respectively. 
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Figure  77:  External  Pressure  Coefficients  for  selected  HALE  Pareto  Members  -  Upper  Surface 
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Figure  78:  External  Pressure  Coefficients  for  selected  HALE  Pareto  Members  -  Lower  Surface 


The  pressure  coefficients  calculated  by  PANAIR  about  the  Pareto  wings  are  of  a  similar  order  of  magnitude 
as  those  computed  for  the  benchmark  wing. 

Figures  79  to  81  show  the  deflected  Pareto  wings.  The  original  wing  position  and  internal  structure  is 
shown  as  a  wire  frame  for  comparison. 
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MSC  Patran  2005  r2  1 2>Jurv06  22  1 9  04 

Fringe  DEFAULT  SCI  A1  Non-linear  100  %  of  Load.  Displacements.  Translational  Magnitude.  (NON-LAYERED) 
Deform  DEFAULT  SCI.  A1  Non-linear  100  %  of  Load.  Displacements.  Translational.  26*000 


Figure  79:  HALE  Pareto  0  Displacement 
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Figure  80:  HALE  Pareto  7  Displacement 
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MSC.Patran  2006  r2  12-Jui>06  22:33:33 
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Figure  81:  HALE  Pareto  23  Displacement 


From  the  above  Figures  showing  the  displacement  of  the  HALE  UAV  wings  it  is  quite  evident  that  the 
wings  deflected  more  than  the  twenty  percent  span  value  indicated  in  Table  15.  Further  testing  is  needed  to 
analyse  the  effects  of  constraining  the  displacement. 


An  analysis  of  the  strains  produced  by  the  aerodynamic  forces  is  shown  in  Figures  82  to  84. 
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Figure  82:  HALE  Pareto  0  Strain 
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Figure  83:  HALE  Pareto  7  Strain 
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Figure  84:  HALE  Pareto  14  Strain 
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Due  to  the  large  displacements  shown  in  Figures  82  to  84,  the  internal  strains  in  the  HALE  UAV  wings  are 
very  large  with  Pareto  0  failing  the  maximum  strain  allowable  value  specified  in  Table  14. 

The  results  for  the  HALE  UAV  wing  optimisation  are  summarised  in  Table  1 6  for  comparison. 
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Pareto  0 

Pareto  7 

Pareto  14 

Benchmark 

Deflection  (%) 

29.65 

31.08 

37.16 

16.7 

Maximum  Strain 

0.00691 

0.000893 

0.00106 

0.000487 

Mass  (kg) 

963.38 

774.28 

525.51 

1138.15 

L/D 

40.805 

36.559 

36.197 

31.729 

Mass  Penalty  (kg) 

286.55 

235.65 

167.69 

0 

Drag  Penalty 

0.0051 

0.0074 

0.0076 

0 

Table  16:  Summary  of  Displacement,  Strain  and  Mass  for  Selected  HALE  Pareto  members 


All  the  Pareto  members  selected  from  the  Pareto  front  shown  in  Figure  68  had  better  Lift  to  Drag  ratios. 

The  mass  of  the  simulated  Pareto  members  was  lower  than  that  of  the  benchmark,  but  the  deflection  of  the 
wings  was  larger  than  the  allowable  twenty  percent.  The  maximum  strain  values  calculated  by 
MSC.Nastran"  were  higher  than  those  calculated  for  the  MALE  UAV  Pareto  wings.  This  is  caused  by  the 
large  displacements  of  the  wings  due  to  a  less  rigid  wing  internal  structure. 

Furthermore,  all  the  HALE  Pareto  wings  had  penalties.  This  would  have  influenced  the  HAPMOEA  code 
as  it  determined  new  generations  of  candidate  solutions  for  evaluation. 

Pareto  0  failed  the  maximum  strain  criteria  as  the  maximum  strain  value  was  almost  twice  the  allowable 
amount.  This  indicates  that  panels  within  the  wing  would  fail  due  to  the  aerodynamic  loads  placed  on  the 
wing  skin. 
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8  Conclusions 


This  report  has  described  the  characteristics  and  implementation  details  of  a  framework  in  which  different 
aeronautical  problems  can  be  analysed.  The  framework  is  based  on  a  Multi-objective  Evolutionary 
Algorithm.  The  report  indicates  and  addressed  some  of  the  limitation  of  current  optimisation  techniques 
and  the  need  for  alternative  approaches.  The  report  provides  a  description  of  the  main  components  of  the 
framework,  i.e.  an  Aero-structural  solver,  a  wing  mesh  generation  capability,  a  robust  MOEA  optimiser  and 
capabilities  for  parallel  computing  and  post  processing.  Hence  we  have  within  the  framework,  a  complete 
set  of  numerical  tools  for  handling  mathematical  test  functions  and  real  world  problems  in  aeronautics. 

The  methodology  is  illustrated  on  its  application  to  the  design  optimisation  of  UAV  wings  as  a  multi¬ 
objective  and  MDO  problem  showing  the  benefits  of  the  method.  The  method  is  found  to  be  capable  of 
identifying  the  trade-off  between  the  multi-physics  involved  and  provides  classical  aerodynamic  shapes  as 
well  as  alternative  configurations  from  which  the  designer  can  choose.  It  is  observed  that  there  was  a 
computational  gain  on  using  asynchronous  evaluation  and  a  hierarchical  topology  of  fidelity  models  as 
compared  to  a  single  model  during  the  optimisation. 

Large  gains  in  aerodynamics  can  be  made  for  both  the  MALE  and  HALE  UAV  test  cases  through  the  use 
of  the  data  calculated  by  the  HAPMOEA  code  from  the  assumed  benchmark  wings.  The  lift  to  drag  ratio 
can  be  extended  to  regions  only  seen  previously  within  the  flight  envelopes  of  high  performance  gliders. 

Structural  gains  can  also  be  made,  though  as  the  true  benchmark  structures  are  more  complex  and  probably 
include  hard  points  within  the  wings  and  fuel  tanks,  the  results  obtained  through  the  HAPMOEA  code  are 
only  indicative  of  future  possibilities  given  more  accurate  data. 

A  more  detailed  analysis  needs  to  be  performed  on  algorithms  to  better  describe  the  internal  structure  of 
high  aspect  wings.  This  analysis  could  aid  the  optimisation  of  the  HALE  UAV  test  case  and  limit  the 
penalties  due  to  deflection  and  strain. 

Although  the  implemented  penalty  functions  and  values  worked  when  optimising  both  the  MALE  and 
HALE  UAV  wings,  a  more  thorough  investigation  of  the  effect  of  different  penalty  functions  and  the 
calculation  of  penalty  values  on  the  optimiser  simulation  solution  could  yield  better  algorithms  which, 
when  used  in  the  HAPMOEA  code,  could  result  in  quicker  solution  times  with  results  which  didn’t  fail  any 
of  the  constraints. 

Given  more  accurate  benchmark  information,  a  more  complex  and  detailed  investigation  and  optimisation 
could  be  performed. 
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9  Future  Work 


A  limitation  to  this  study  is  that  the  solution  found  is  not  a  full  aero-elastic  solution  but  rather  a  one 
iteration  coupling  of  the  aerodynamics  of  the  candidate  wing  with  the  structural  model.  Furthermore,  a 
rather  low  fidelity  model  is  used  to  calculate  the  form  and  skin  friction  components  of  the  drag  produced  by 
the  candidate  wing. 

Further  analysis  methods  and  techniques  could  be  incorporated  into  the  Aero- Structural  Solver  framework, 
or  modifications  to  the  code  could  be  made  to  perform  calculations  which  would  assist  in  designing  more 
optimal  wings  based  upon  other  constraints.  The  work  could  be  divided  into  two  sections  namely  the 
aerodynamics  and  aerodynamic  package,  and  the  structural  model  and  structural  simulation  method. 

9.1  Aerodynamics 

9.1.1  Solver 

The  PANAIR  aerodynamic  solver  could  be  upgraded  to  make  use  of  a  Navier-Stokes  solver.  This  would 
increase  the  computational  cost  for  each  candidate  wing  simulation,  but  would  also  increase  the  accuracy  of 
the  solution  attained  if  complex  flow  geometries  were  being  analysed  such  as  slotted  flaps. 

If  the  solver  were  upgraded,  special  consideration  would  have  to  be  given  when  designing  the  program  to 
automatically  generate  the  domain  meshes  about  the  planned  wings.  If  little  consideration  was  given  in 
defining  mesh  densities  about  areas  of  high  macroscopic  gradients,  loss  of  accuracy  could  be  significant. 

9.1.2  Drag  Calculations 

Currently  the  Aero- Structural  Solver  makes  use  of  PANAIR  and  FRICTION  to  calculate  the  aerodynamic 
loads  and  aerodynamic  fitness  values  for  each  candidate  wing.  As  PANAIR  is  only  applicable  to  linear 
potential  flows,  it  would  be  beneficial  to  upgrade  to  a  nonlinear  flow  solver  such  as  TRANAIR®.  Such  an 
upgrade  would  allow  the  user  or  optimiser  the  ability  to  model  complex  flows  which  included  shock  waves 
at  or  near  transonic  speeds.  This  would  also  minimise  any  solution  errors  due  to  large  Mach  angles. 

FRICTION  is  a  low  fidelity  program  which  was  passed  a  constant  value  for  the  percentage  chord  location 
where  the  transition  from  a  laminar  to  turbulent  boundary  layer  took  place.  The  assumption  was  a  crude 
estimate  made  by  the  authors  and  would  not  hold  true  for  all  the  candidate  wings  analysed.  For  a  more 
accurate  skin  friction  drag  coefficient  calculation,  a  higher  fidelity  model  could  be  incorporated  into  the 
solution  sequence.  Another  simple  addition  which  could  be  made  to  FRICTION  could  be  a  simple  iterative 
method  which  could  determine  the  Reynolds  number  at  which  the  flow  would  make  the  transition  from 
laminar  to  turbulent,  and  the  corresponding  chord  ratio  passed  to  FRICTION  when  calculating  the  required 
coefficients  of  drag. 

9.2  Structures 

9.2.1  Materials 

MSC.Nastran8  has  the  ability  to  define  materials  that  are  made  up  of  a  number  of  plies.  Hence  graphite 
epoxy  based  materials  could  be  correctly  defined  through  the  build  up  and  orientation  of  the  plies  within 
the  material.  This  could  lead  to  a  more  detailed  optimization  where  the  material  thickness  could  not  only  be 
a  variable,  but  the  make  up  of  the  material,  fibre  thickness,  orientation  and  laminate  layers,  could  also  be 
added  as  design  variables  in  the  optimization  process. 
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9.2.2  Elements 


Most  of  the  element  input  decks  used  when  defining  the  elements  in  the  structural  model  had  the  capacity 
to  make  use  of  material  orientation  within  the  element.  Future  work  could  entail  the  use  of  fibre  orientation 
material  properties  such  as  those  mentioned  above.  This  would  lead  to  a  better  understanding  and 
optimisation  as  the  element’s  true  behaviour  under  loading  could  be  studied. 

If  the  position  of  the  spars  and  ribs  as  functions  of  chord  and  span  respectively  were  included  as  design 
variables  within  the  simulation  framework,  the  optimality  of  the  spars  could  be  increased  as  spars  could  be 
placed  such  that  the  bending  moments  generated  by  the  aerodynamic  forces  could  be  minimised,  and  the 
variable  rib  placement  could  mean  that  better  rib  distributions  could  be  achieved  where  sections  of  high 
strain  which  occur  at  crank  locations  could  be  minimised  through  the  use  of  two  closely  spaced  ribs. 

More  complex  thickness  functions  could  be  devised,  such  as  third  or  fourth  order  polynomials,  to  describe 
the  way  in  which  Spar  and  Rib  thicknesses  vary  according  to  span.  This  could  assist  in  the  reduction  of 
high  strain  concentrations  which  would  occur  where  the  wing  geometry  changes  greatly. 

CBAR  or  CBEAM  elements  could  be  used  to  describe  the  Rib  and  Spar  Caps.  The  height  and  width  of  the 
sections  could  then  be  added  to  the  design  variables  optimised  by  the  Evolutionary  Optimiser  code.  This 
would  allow  for  more  accurate  deflection  calculations  and  assist  in  the  minimisation  of  the  candidate  wing 
structure. 

9.2.3  Solution  Method 

Further  analysis  could  be  conducted  to  incorporate  buckling  analysis  to  investigate  the  natural  frequencies 
of  the  candidate  wings.  MSC.Nastran8  can  incorporate  the  MSC. Flightloads “  package  and  undertake  full 
aero-elastic  analysis  of  candidate  wing  geometries.  This  would  greatly  increase  the  solution  resolution  of 
the  candidate  wings  as  solutions  could  be  penalised  against  a  final  wing  deflection  value  instead  of  the  first 
solution  iteration.  This  is  beneficial  as  many  wing  solutions  would  solve  to  a  position  which  would  lie  in 
the  predefined  constraints,  but  as  the  initial  iteration  indicates  that  the  deflection  is  greater  than  this 
predefined  value,  the  wing  is  discarded  prematurely. 
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